Erwan Drezen

Guillaume Rizk

Rayan Chikhi

Charles Deltel

Claire Lemaitre

Pierre Peterlongo

Dominique Lavenier

GATB: Genome Assembly & Analysis Tool Box

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The analysis of NGS data remains a time and space-consuming task. Many efforts have been made to provide efficient data structures for indexing the terabytes of data generated by the fast sequencing machines (Suffix Array, Burrows-Wheeler transform, Bloom Filter, etc.). Genome assemblers such as Velvet (Zerbino et al., 2008), ABySS (Simpson et al., 2008), SOAPdenovo2 (Luo et al., 2012), SPAdes (Bankevich et al., 2012) or mappers such as BWA (Li et al., 2009) or variant detection such as CRAC (Philippe et al., 2013) for instance make an intensive use of these data structures to keep their memory footprint as low as possible.

At the same time, parallelism has been largely investigated to reduce execution time. Many strategies such as GPU implementation (Liu et al., 2012), cloud deployment (Zhao et al., 2013), algorithm vectorization (Rizk et al., 2010), multithreading, etc., have demonstrated high potentiality on NGS processing.

The overall efficiency of NGS software depends on a smart combination of data representation and use of the available processing units. Developing such software is thus a real challenge, as it requires a large spectrum of competence from high-level data structure and algorithm concepts to tiny details of implementation.

The GATB library aims to ease the design of NGS algorithms. It offers a panel of high-level optimized building blocks to speed-up the development of NGS tools related to genome assembly and/or genome analysis. The underlying data structure is a memory efficient de-Bruijn graph (Compeau et al., 2011), and the general parallelism model is multithreading. The GATB library targets standard computing resources such as current multicore processor (laptop computer, small server) with a few GBytes of memory.

Hence, from the high-level C++ functions available in the GATB library, NGS programing designers can rapidly elaborate their own software based on state-of-the-art algorithms and data structures of the domain.

Based on the same idea, other bioinformatics libraries exist, from which domain specific tools can be elaborated. The NGS++ library (Markovis et al., 2013) is specifically tailored for developing applications that work with genomic regions and features, such as epigenomics marks, gene features and data that are associated with BED type files. The SeqAn library (Doring et al., 2008) is a general-purpose library targeting standard sequence processing. Advanced data structures such as de-Bruijn graphs are not included in SeqAn. Khmer (Crusoe et al., 2014) is a library and toolkit for doing k-mer-based NGS dataset analysis. As with GATB, most of khmer relies on an underlying probabilistic data structure (Bloom filter). The khmer library can be used in various k-mer processing such as abundance filtering, error trimming, graph size filtering or partitioning.

Methods

One of the main concerns of the GATB-core library is to provide computing modules able to run on standard machines, i.e. computers not requiring large amount of main memory. The central data structure is a de-Bruijn graph from which numerous actions can be performed: data error correction, assembly, biological motif detection (e.g. SNP), etc. The graph is constructed by extracting and by counting all the different k-mers from one or several sequencing data sets. This time and space consuming task is conducted by a disk streaming algorithm, DSK, (Risk et al., 2013) which adapts its memory requirement according to the available computer memory. Tradeoff between execution time and memory occupancy can be set up: the larger the computer memory, shorter the computation time (reduced disk access).

The de-Bruijn graph memory footprint is kept very low thanks to an optimized Bloom filter representation (Chikhi et al., 2012;Salikhov et al., 2013). Only vertices of the de-Bruijn graph are memorized. Edges are deduced by querying the Bloom filter. False positives (due to the probabilistic behavior of the Bloom filter) are suppressed by adding an extra data structure enumerating critical vertices. This very efficient de-Bruijn graph representation fits, for example, a complete mammal genome in about 4 GB.

Implementation

The GATB library is composed of five main packages: system, tools, bank, kmer and de-Bruijn packages. The system package holds all the operations related to the operating system: file management, memory management and thread management. Using such an abstraction allows client code to be independent from the OS, thus suppressing compilation directive inside the code or improving some OS accesses by hiding specific OS optimization. The supported operating systems are Linux, Mac and Windows.

The tools package offers generic operations used throughout user application, but not specific to genomics area. For example, this package includes design pattern tools (such as iterators, observers, smart pointers, etc.) and object collections (such as containers, bags, iterables, etc.). It also optimizes the way GATB data structures are saved. The HDF5 file format is currently used (HDF5 2012). This powerful technology is extremely well suited for large and complex data collection such as those handled in the GATB library. The bank package provides operations related to standard genomic sequence dataset management. All the main sequence file formats are supported, and high-level interfaces allow sequences to be easily iterated regardless of the input format. In other words, algorithms are written independently of the input formats. The kmer package is dedicated to fine-grained manipulation of k-mers. Optimized routines are provided to perform k-mer counting from large sequence datasets, to find k-mer neighborhood or to select k-mers based on different criteria.

Finally, the de-bruijn package provides high-level functions to manipulate a static de-Bruijn graph data structure: creation from a set of k-mers, iteration through different nature of nodes (simple k-mers, branching k-mers, etc.), extraction of neighbor nodes, etc. Additional information (e.g. k-mer coverage, markers of visited nodes) is stored in the graph branching nodes. From this abstraction level, developing new tools based on de-Buijn graphs is fast, and does not require programmers to delve into low level details.

The GATB library takes benefit of the parallel nature of today's multicore architecture of microprocessors. When possible, time-consuming parts of the code are multithreaded to provide fast runtime execution. The GATB library is developed in C++ under the A-GPL license and is available from the following web site: http://gatb.inria.fr. An extensive documentation with tutorials is available to guide designers in the process of developing new NGS tools from the GATB building blocks: http://gatb-core.gforge.inria.fr (see also supplementary file 2 for technical implementation details).

Result

To demonstrate the efficiency of the GATB library, a few software implemented from GATB are briefly presented. The idea is to give a quick overview of the application spectrum of the GATB library, and some performance numbers.

Minia (Chikhi et al., 2012) is a short-read de-Bruijn assembler capable of assembling large and complex genomes into contigs on a desktop computer. The assembler produces contigs of similar length and accuracy to other de-Bruijn assemblers -e.g. Velvet (Zerbino et al., 2008). As an example, a Boa constrictor constrictor (1.6 Gbp) dataset (Illumina 2x120 bp reads, 125x coverage) from Assemblathon 2 (Bradman et al., 2013) can be processed in approximately 45 hours and 3GB of memory on a standard computer (3.4 GHz 8-core processor) using a single core, yielding a contig N50 of 3.6 Kbp (prior to scaffolding and gap-filling).

Bloocoo is a k-mer spectrum-based read error corrector, designed to correct large datasets with very low memory footprints. It uses the disk streaming k-mer counting algorithm contained in the GATB library, and inserts solid k-mers in a bloom-filter. The correction procedure is similar to the Musket multistage approach (Liu et al., 2012). Bloocoo yields similar results while requiring far less memory: for example, it can correct whole human genome re-sequencing reads at 70 x coverage with less than 4GB of memory (see supplementary file 1 for extra information on Bloocoo).

DiscoSNP aims to discover Single Nucleotide Polymorphism (SNP) from non-assembled reads and without a reference genome. From one or several datasets a global de-Bruijn graph is constructed, then scanned to locate specific SNP graph patterns (Uricaru et al., 2014). A coverage analysis on these particular locations can finally be performed to validate and assign scores to detected biological elements. Applied on a mouse dataset (2.88 Gbp, 100 bp Illumina reads), DiscoSnp takes 34 hours and requires 4.5 GB RAM. In the same spirit, the TakeABreak software discovers inversion variants from non-assembled reads. It directly finds particular patterns in the de-Bruijn Graph, and provides execution performances similar to DiscoSNP (Lemaitre et al., 2014).

• Liu Y, Schröder, J. and Schmidt, B ( 2013

Figure 1 :

 1 Figure 1: Schematic view of the GATB organization

  ) Musket: a multistage k-mer spectrum-based error corrector for Illumina sequence data, Bioinformatics, 29, 308-315 • Liu Y, Schmidt B, Maskell D (2012) CUSHAW: a CUDA compatible short read aligner to large genomes based on the Burrows-Wheeler transform. Bioinformatics Volume 28, Issue 14, Pp. 1830-1837 • Luo R, Liu B, Xie Y, Li Z, Huang W, et al. (2012) SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience 1: 18 • Markovits et al. (2013) NGS++: a library for rapid prototyping of epigenomics software tools. Bioinformatics, 29 (15): 1893-1894 • Philippe N, Salson M, Commes T, Rivals E (2013) CRAC: an integrated approach to the analysis of RNAseq reads. Genome Biology 14:R30 • Rizk G, Lavenier D (2010) GASSST: global alignment short sequence search tool. Bioinformatics, vol 26, issue 20, pp 2534-2540 • Rizk G, Lavenier D, Chikhi R (2013) DSK: k-mer counting with very low memory usage. Bioinformatics 29 (5), 652-653 • Salikhov K, Sacomoto G, Kucherov G (2013) Using cascading bloom filters to improve the memory usage for de-Bruijn graph, Algorithms in Bioinformatics, Lecture Notes in Computer Science Volume 8126, pp 364-376 • Simpson JT, Wong K, Jackman SD, Schein JE, Jones SJM, et al. (2009) ABySS: a parallel assembler for short read sequence data. Genome Res 19: 1117-1123 • Uricaru R, Rizk G, Lacroix V, Quillery E, Plantard O, Chikhi R, Lemaitre C, Peterlongo P (2014) Referenceree detection of genotypable SNPs, in revision to NAR • Zhao S and al. (2013) Rainbow: a tool for large-scale whole-genome sequencing data analysis using cloud computing. BMC Genomics 2013, 14:425 • Zerbino, DR., Birney E (2008) Velvet: algorithms for de novo short read assembly using de-Bruijn graphs. Genome research 18.5: 821-829

Acknowledgment

Funding: this work is supported by the ANR (French National Research Agency), ANR-12-EMMA-0019-01

• Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, et al. (2012) SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol 19: 455-477.


Erwan Drezen
email: erwan.drezen@inria.fr

Dominique Lavenier

Quality metrics for benchmarking sequences comparison tools

Keywords: sequence, alignment, benchmark

Background: Comparing sequences is a daily task in bioinformatics and many software try to fulfill this need by proposing fast execution times and accurate results. Introducing a new software in this field requires to compare it to recognized tools with the help of well defined metrics.

Results: A set of quality metrics is proposed that enables a systematic approach for comparing alignment tools. These metrics have been implemented in a dedicated software, allowing to produce textual and graphical benchmark artifacts.

Introduction

In Bioinformatics, the task of comparing genomic sequences is ubiquitous, leading to various software proposals over the years. With the NGS revolution, this task is becoming more and more time consuming as the size of NGS data truly explodes.

Software have tried to cope with this situation in different ways. First, by using known algorithms with optimized hardware usage. Second, by proposing new heuristics to reduce the searching space. The historically sequence comparison tools are:

• SSEARCH [START_REF] Pearson | Improved tools for biological sequence comparison[END_REF] that proposes an implementation of the Smith-Waterman algorithm with an efficient usage of multicore / SSE architecture; it provides an "exact" search alignment algorithm and is therefore slow, and consequently inadapted for large requests. • BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF] and FASTA [START_REF] Pearson | Improved tools for biological sequence comparison[END_REF] that both propose a seed based heuristics to speedup the computation with a small loss of quality w.r.t. exact algorithms.

The trade-off between speed and quality is the key point. As a consequence, to evaluate new Alignement Search Tools (AST for short), benchmarking has to focus on time and quality metrics. Comparing execution times is straightforward. On the other hand, comparing quality is more challenging since it requires to define a precise quality metric.

The paper presents a methodology for comparing the results of AST. We explicitely target intensive sequence comparison for which there is an increasing demand due to high throuhput sequencing oportunities. A quality metric is defined and applied to a set of software and benchmarks.

Benchmarking AST

Configuration

To fairly benchmark AST, many configuration aspects must be taken into consideration:

1 For a given request, AST have to be configured for providing similars results. This is far from obvious since they generally propose a large set of parameters.

In some cases, they can have hidden configuration. 2 Some AST implicitely modified the input banks. For instance, low complexity region [START_REF] Wootton | Statistics of local complexity in amino acid sequences and sequence databases[END_REF], [START_REF] Morgulis | Fast and Symmetric DUST Implementation to Mask Low-Complexity DNA Sequences 5[END_REF] can be systematically masked. These transformations may have significant influence on the results compared to other tools that are not applying these kinds of pre-processing. 3 AST often provide a statistical model telling whether an alignment has to be kept or not (for instance Karlin/Altschul model of Blast). Thus, AST with different models are likely to produce different alignments.

These points have a strong impact on the alignment list generated by the AST. Consequently, they need to be included in the quality metric itself. As a side effect, the AST command lines have to be explicitely described as part of the benchmark results.

In the following, we note T (0) (or T for short) an AST with its default configuration. Specific configurations of T will be noted as T (i) for i ≥ 1.

Input request

To compare AST, a minimal set of common parameters needs to be specified. The following table lists the parameters we selected: Parameters r 1 to r 9 have a direct impact on the output alignments. The r 10 parameter has an impact on the execution time. The execution time can be used to analyze the scalability of AST as a function of the number of cores.

A request R can now be defined by a set of r i . If not all r i are specified, it is supposed that default values are used instead. A relation of equivalence can also be defined between two requests x and y:

x ∼ y ⇐⇒ ∀i ∈ [1..10]; r i (x) = r i (y)
Here is an example of a request R 1 :

R 1 =      r 1 = "swissprot" r 2 = "yeats" r 3 = 1e -5 2.

Output data

Comparing AST qualities implies that AST generates the same nature of information. We choose the information provided by the Blast tabular output format; for each alignment, it provides 12 properties p i as shown in table 2. Table 3 is an illustration of a Blast output. Relying only on this set of properties means that only AST able to provide this information would be candidate. Fortunately, many AST can be parametrized to generate Blast-like tabular output format. 3 Alignment Quality Metric

Definitions

Before specifying an alignment quality metric, we introduce the following definitions:

• I is a set of integers i ∈ [1.
.12]. I specifies a subset of properties p i defined in table 2 • For a given set I, an alignment a is an object having properties p i for i ∈ I • Two alignments x and y are equivalent if their properties have the same values:

x ∼ y ⇐⇒ ∀i ∈ I; p i (x) = p i (y)
• an alignment set A contains alignments defined for a specific set I An alignment set may contain equivalent items. In such a case, items equivalent to some alignment a are grouped into an equivalence class ã.

ã = {x ∈ A; a ∼ x}
We define à as the set built from A by keeping only one representant of each equivalence class ã.

Taking the example of table 3, for I = {1}, we have the alignment set A = {a1, a2, a3} and the set of equivalent alignments à = { ã1 , ã2 } with: With these definitions, an AST output can be formally transformed into an alignment set A where only some properties defined by the set I are kept. An alignment set can thus be defined as A := ϕ(T, R, I) with T the AST, R the request and I the set of properties. Figure 1 illustrates the transformation of a Blast output. Note also that the same output generated by a tool will be viewed differently for two distinct sets I 1 and I 2 . Therefore several quality metrics based on different I sets can be defined to provide different benchmark points of view.

     p 1 (a1) = "SP O13535" p 1 (a2) = "SP O13535" p 1 (a3) = "SP P 02400" p 1 ( ã1 ) = "SP O13535" p 1 ( ã2 ) = "SP P 02400"

Comparing AST

For a request R (defined by a set of r i from table 1) and a given set I, we note:

A 1 = ϕ(T 1 , R, I), alignment set f ound by AST T 1 f or request R A 2 = ϕ(T 2 , R, I), alignment set f ound by AST T 2 f or request R
Comparing T 1 to T 2 consists in finding a set M (I) such as:

M (I) =      (x, y) ∈ A 1 × A 2 ; x ∼ y ∀(x , y ) ∈ A 1 × A 2 ; x ∼ y ⇒ y = y x ∼ y ⇒ x = x
The intuitive idea of this definition is to map an item from one set to at most one item of the other set. Note that some alignments in one set may not be connected to an alignment in the other set as shown figure 2.

The estimation of "how close two sets A 1 and A 2 are" can be done by considering the cardinal |M (I) | which represents the number of common alignments between A 1 and A 2 . We define three numbers N Finally, the alignement quality metric for a given I by the Q (I) function is defined as:

A 1 × A 2 Q (I) ---→ N (I) c , N (I) 1 , N (I) 2
With the following properties:

1

Q (I) (A, A) = (|A|, 0, 0) 2 Q (I) (A 1 , A 2 ) = (x, y, z) =⇒ Q (I) (A 2 , A 1 ) = (x, z, y)
Note that this definition of the alignement quality metric is relative: it only gives information for a specific couple of alignment sets. It does not provide absolute assessment for one AST. The consequence is that one of the AST (T 1 , T 2 ) needs to be chosen as a reference when comparing many AST together.

The percentage of alignments found by T 2 in the alignments of T 1 can also be defined as:

α (I) 2 = N (I) c N (I) c + N (I) 1
Taking figure 2 as example, we have Q (I) (A 1 , A 2 ) = 5, 11, 6 . T 2 is able to find α Comparing N AST T i is simply done by computing Q (I) (A i , A j ) for all (i, j) ∈ [1..N ] 2 . Note that computation is only required for i >= j because of the second property of the quality definition.

Results

List of benchmarked tools and quality metrics

The evaluation of our methodology has been done with the following tools :

T 1 :=SSEARCH T 2 :=BLAST T 3 :=PLAST T 4 :=UBLAST
In this set of AST, SSEARCH is the only exact tool, the other ones are based on heuristics. Thus, SSEARCH is taken as the reference. In that way, it can be seen how close heuristics based tools are from an exact reference.

Heuristics based tools are configured by default to get the best trade-offs between speed and quality. Of course, these tools can also be configured by setting specific parameters in order to get a better quality or a better speed. Table 4 gives several configurations used for the benchmark. Three different quality metrics are considered:

• I query = {1} with Ã, the equivalent class alignment set; this metric compares the number of matched queries between tools • I hit = {1, 2} with Ã, the equivalent class alignment set; this metric compares the number of matched couples [subject,query] between tools • I align = {1, 2, 7, 8, 9, 10} with A, the alignment set; this metric compares the number of alignments sharing the same [subject,query] identifiers and the same [subject,query] boundaries. These 3 metrics go from the worst to the best precision. As a matter of fact, the correponding sets I provide more and more selective relations of equivalence between alignments. As explained later, two AST may be close for I query and even for I hit but may have significant differences for I align ; this is shown by small value for N

Results

In the following tables, we dump α (I) 2 , the percentage of alignments found by T 2 in the alignment set of T 1 .

Two ratios of execution time between T 1 and T 2 can be defined: 1 SU total : speedup for the total execution time of a request, including bank preparation when needed (use of makeblastdb for instance) 2 SU AST : speedup for the execution time of the sequences comparison only (exclude makeblastdb execution time for instance) Tables 5 and6 show a protein/protein comparison between swissprot [5] and two bacterias [START_REF]Proteome of Escherichia Coli[END_REF], [START_REF]Proteome of Chitinophaga Pinensis[END_REF] from which several comments can be done:

• Heuristics based tools have good speedups, especially UBLAST • AST quality are very close to SSEARCH for I (query) (except UBLAST, Table 6) • With default configuration, BLAST and UBLAST have poor results for I (hit) and I (align) ; only PLAST manages to recover most of the SSEARCH results • For two similar requests, PLAST has the smallest discrenpencies for the 3 metrics, followed by BLAST and then by UBLAST Tables 7 compares a subset of 40000 sequences from swissprot to another subset of 40000 sequences from swissprot. Here, BLAST and PLAST have similar values for the 3 metrics and are close to the reference; speedups are interesting, especially for PLAST. For this request, UBLAST is still very fast but has poor values for the 3 metrics. Other tools configurations give similar quality results. Tables 8 shows a genomic query and table 9 a meta-genomic query [START_REF]Environmental Samples from TARA[END_REF], with swissprot as subjet bank. SSEARCH can't be used as reference here so BLASTQ is used instead. BLAST and PLAST are similar for I (align) , with better results for PLAST on I (hit) and better results for BLAST on I (query) . UBLAST is still the fastest tool with poor quality, in particular in the meta-genomic case. In brief, we summarize this benchmark with the following observations:

• As expected, heuristics based tools are fast at the expense of quality; for instance, UBLAST is very fast but has often poor quality for the chosen metrics

• Surprisingly, BLAST with default configuration has sometimes poor quality for some metrics; a special configuration like BLASTQ is needed for recovering full quality • PLAST is a good trade-off between speed and quality; it is faster than BLAST with similar quality values.

Conclusion

In this paper, we have defined a set of quality metrics for comparing alignment search tools (AST). Comparing AST is important to understand how the AST behave and what should be expected from the alignments they produce. In particular, these metrics allow the distance between exact and heuristics based tool to be precisely evaluated. These metrics have been tested on standard AST. The benchmark results provides some interesting hints that can be used to match a specific tool with user needs; for instance, PLAST is a good trade-off in terms of speed/quality, and UBLAST is interesting if speed is crucial and if quality loss is acceptable (particularly for alignment and hit metrics).

Other quality metrics are currently under investigation to better capture user needs according to the application domains where AST are involved. A next step is also to gather in a public database many experimentatons and offer to the scientific community a way to select the best AST according to their needs. Other AST will also be added to extend the scope of this methodology. Finally, a web interface [START_REF]Web Interface for Metrics Display[END_REF] is under development and already provides a graphical representation of the metrics.

  c := |M (I) |, number of common alignments between A 1 and A 2 • N (I) 1 := |A 1 | -|M (I) |, number of alignments specific to A 1 • N (I) 2 := |A 2 | -|M (I) |, number of alignments specific to A 2

Figure 2

 2 Figure 2 Mapping equivalent alignments between two sets A 2 A 1

  .25% of the alignments found by T 1 . Similary, T 1 finds α .45% of the alignments found by T 2 .

Table 1

 1

	Request parameters
	r 1 ) query bank	r 6 ) reward cost
	r 2 ) subject bank	r 7 ) penalty cost
	r 3 ) e-value	r 8 ) substitution matrix
	r 4 ) open gap cost	r 9 ) low complexity filtering
	r 5 )

extend gap cost r 10 ) number of cores

Table 2

 2 ) percentage of identity p 9 ) start offset in subject p 4 ) alignment length p 10 end offset in subject p 5 ) number of mismatches p 11 ) e-value p 6 ) number of gap openings p 12 ) bitscore

	Alignment properties
	p 1 ) query sequence id	p 7 ) start offset in query
	p 2 ) subject sequence id	p 8 ) end offset in query
	p 3

Table 3

 3 Example of Blast tabular output

	p 1	p 2	p 3	p 4	p 5	p 6	p 7	p 8	p 9	p 10	p 11	p 12
	SPO13535 SPP02309	39.34 305 181 4	29 332 17 302 7.5e-62 231.6
	SPO13535 SPQ91G55	34.93 209 134 2	42 250 46 252 1.5e-31 130.4
	SPP02400	SPQ196T6 30.18 603 388 3	90 663 49 647 1.6e-76 282.4

  Figure 1 Building an alignment set from an AST output

	r2 : yeats Request r1 : swissprot	BLAST	Alignments Properties I = {1, 2, 7, 8, 9, 10}
	r3 : 1e-5
		BLAST Output	Alignments Set
	SPO13535	SPP02309	39.34	...
	SPO13535	SPQ91G55	34.93	...
	........	........	.....	...
	3.2 From AST output to alignment set

Table 4

 4 Tools specific configurations

	Name	Configuration	Purpose
	BLASTN	blast -task blastn	nucl/nucl requests
	MEGABLAST	blast -task megablast	nucl/nucl requests
	BLASTQ	blast -max-target-seqs 25000 quality improvement
	UBLASTQ	ublast -accel 1	quality improvement
	PLASTS	plast -seeds-use-ratio 0.01	speed improvement

Table 5 [

 5 Escherichia Coli vs. uniprot sprot, evalue=1e-3, nbcores=16]

	T 1	T 2	α 2 (I align ) α 2 (I hit )	α 2 (Iquery)	SU total	SU AST
	ssearch	blast	55.5	58.1	99.8	5.48	5.80
	ssearch	plast	93.9	97.1	99.8	26.22	32.54
	ssearch	ublast	29.0	64.4	99.8	64.41	146.60
	ssearch	blastQ	94.2	97.4	99.8	4.90	5.08
	ssearch	plastS	73.2	77.1	99.4	52.07	84.43
	ssearch	ublastQ	30.4	69.1	99.8	40.61	62.63
	Table 6 [Chitinophaga Pinensis vs. uniprot sprot, evalue=1e-3, nbcores=16]
	T 1	T 2	α 2 (I align ) α 2 (I hit ) α 2 (Iquery)	SU total	SU AST
	ssearch blast	61.4	64.2	98.7	5.33	5.45
	ssearch plast	92.0	95.9	97.9	29.23	33.26
	ssearch ublast	19.9	52.4	78.9	89.48	157.96

Table 7

 7

	[uniprot sprot 40000b vs. uniprot sprot 40000, evalue=1e-30, nbcores=16]
	T 1	T 2	α 2 (I align ) α 2 (I hit ) α 2 (Iquery)	SU total	SU AST
	ssearch blast	86.5	89.9	97.1	8.40	8.46
	ssearch plast	88.9	92.4	96.3	44.62	46.48
	ssearch ublast	6.4	19.9	22.8	375.46	591.03

Table 8 [

 8 SRR142736 vs. uniprot sprot, evalue=1e-3, nbcores=16]

	T 1	T 2	α 2 (I align ) α 2 (I hit )	α 2 (Iquery) SU total	SU AST
	blastQ blast	73.8	72.3	100.0	1.00	1.00
	blastQ plast	77.7	95.4	98.1	6.70	7.10
	blastQ ublast	19.8	64.9	82.0	40.30	68.91
	Table 9 [TARA sample vs. uniprot sprot, evalue=1e-3, nbcores=16]
	T 1	T 2	α 2 (I align ) α 2 (I hit )	α 2 (Iquery) SU total	SU AST
	blastQ blast	82.1	80.8	100.0	1.00	1.00
	blastQ plast	80.8	94.7	83.1	6.50	6.68
	blastQ ublast	17.2	59.2	29.9	71.10	105.50

Table 10 [

 10 SRR027344 vs. nt 2000000 seqs, evalue=1e-30, nbcores=16]

	T 1	T 2	α 2 (I align ) α 2 (I hit )	α 2 (Iquery)	SU total	SU AST
	blastQ blast	64.8	64.6	100.0	1.03	1.03
	blastQ megablast	61.1	61.0	96.6	5.17	12.48
	blastQ plast	94.2	94.2	96.8	5.38	14.07
	blastQ ublast	90.1	94.9	95.6	0.63	0.65
	Tables 10 is a nucleotide/nucleotide request, comparing a set of reads to 2, 000, 000
	sequences from the nt [9] database. Here, PLAST and UBLAST provide the best
	quality metrics w.r.t. BLASTQ, with good speedups for PLAST.


Erwan Drezen

Patrick Durand

Dominique Lavenier

Dominique Lavenier Klast
email: klast@korilog.com

KLAST integration

Professional version of PLAST (BMC Bioinformatics, 2009) Optimized for bank-to-bank sequence comparisons

Provide high speed and high quality results

  bacterial genomics KLAST and BLAST benchmark comparison of 8,245 sequences (translated 454 reads) from Tara Oceans metagenomic data against 15 million proteins from Uniprot. Both algorithms ran on 8 Intel Xeon cores. KLAST achieved sequence comparisons 18x times faster than BLAST, while covering up to 96% of the results produced by BLAST. Benchmark data courtesy of Thomas Vannier and Jean-Marc Aury research team (Genoscope/CEA). More on this study is available at tinyurl.com/d54ahrb The SSEARCH, BLAST, USEARCH software are considered Comparison of 2,329 protein sequences from bacterium A.hospitalis against the SwissProt databank The reference is given by SSEARCH since it implements the rigorous Smith and Waterman algorithm and generates optimal alignments. Alignments are evaluated on a moderate-size dataset due to the long execution time of SSEARCH. The diagrams synthetize the numbers of alignments found by the different softwares (A, B), the number of queries matching the reference databank (C) and the search time (D). red section reports other alignments. As it can be seen, the total number of alignments can exceed the number of alignments found by SSEARCH. This is mainly due to the fragmentation of long alignments (including large gaps) into shorter ones by KLAST, BLAST and USEARCH. More benchmarks are available at koriscale.inria.fr Command-line tool Workflow and data analytics platform Graphical platform C L C G e n o m i c s workbench

E

  . DREZEN -P. DURAND -D. LAVENIER Bio--IT -Boston -April 29 th -1 st May 2014


A C Levasseur-Regourd

J Lasue

Comparison between interplanetary and cometary dust, from polarimetric remote and in situ studies

Properties of dust in the interplanetary dust cloud and in cometary comae are approached through polarimetric observations, including polarimetric imaging and analyses of the phase angle and wavelength dependences. Local changes are clues to changes in the properties of the dust particles. They are likely to correspond to modifications in the size distribution and to evolution processes related to the evaporation (or alteration) of ices and organics partly constituting some dust particles.

Intensity and linear polarization observations

Linear polarization of solar light scattered by dust particles forming clouds or regolitic surfaces depends on the phase angle and the wavelength of the observations, as well as on the intrinsic properties of the scattering medium. In the visible domain, phase curves of comets and zodiacal light suggest that the particles are irregular, with sizes mostly above one micron. Intensity and polarization observations provide information on dust properties, such as spatial density, local changes, size distribution, morphology, composition and albedo. Observations of solar light scattered by interplanetary dust, so-called the zodiacal light, have been used to derive tables and maps in helio-ecliptic coordinates, with a resolution of 5° or larger [START_REF] Levasseur-Regourd | Absolute photometry of zodiacal light[END_REF][START_REF] Leinert | The 1997 reference of diffuse night sky brightness[END_REF][START_REF] Levasseur-Regourd | Optical and thermal properties of interplanetary dust[END_REF][START_REF] Lasue | Interplanetary dust. Polarimetry of stars and planetary systems[END_REF]. They are critical for studies of extended faint astronomical sources in the background [START_REF] Levasseur-Regourd | Optical and thermal properties of interplanetary dust[END_REF]. Besides, they lead to local information from partial inversion of the line of sight integrated data [START_REF] Leinert | The 1997 reference of diffuse night sky brightness[END_REF]. The local polarization at 90° phase angle increases with increasing solar distance, at least up to 1.5 au in the ecliptic plane, while the local albedo decreases [START_REF] Leinert | The 1997 reference of diffuse night sky brightness[END_REF][START_REF] Lasue | Interplanetary dust. Polarimetry of stars and planetary systems[END_REF]. In cometary comae, evidence for changes in local polarization with the distance from the nucleus and with wavelength (through narrow continuum filters) had been discovered through observations of the OPE instrument on-board the Giotto spacecraft during its flyby of comet 1P/Halley [START_REF] Levasseur-Regourd | Similarities between in situ measurements of local dust scattering and dust flux impact data within the coma of 1P/Halley[END_REF]. Since then, numerous remote observations in the visible and near infrared domains have provided clues, through polarimetric imaging techniques, to heterogeneities in the dust properties within cometary comae, with an increase in jet-like polarization signatures (as noticed for 67P/Churyumov-Gerasimenko, the target of Rosetta mission, after its 2009 perihelion passage) and a possible decrease of the polarization in the nearnucleus region [e.g., [START_REF] Hadamcik | Light scattering by irregular dust particles in the solar system: observations and interpretation by laboratory measurements[END_REF][START_REF] Kolokolova | Physical properties of cometary dust from light scattering and thermal emission[END_REF][START_REF] Hadamcik | Polarimetric observations of comet 67P/Churyumov-Gerasimenko during its 2008-2009 apparition[END_REF]. It was noticed that polarization phase curves provide clues to different type of dust particles in cometary comae, and that the polarization at a given phase angle usually increase with increasing wavelength, at least in the visible domain.

Interpretation of observations

Analysis of such studies are revisited, both for interplanetary and cometary dust, with emphasis on the significance of complementary experimental and numerical simulations to provide relevant interpretations. Interplanetary dust and cometary dust particles are noticed to present some significant similarities [START_REF] Levasseur-Regourd | Physical properties of cometary and interplanetary dust[END_REF]. Numerical simulations of the polarimetric behaviour of interplanetary dust particles suggest that, in the inner solar system, they consist of absorbing (organic compounds) and less absorbing (silicates) materials, that radial changes originate in a decrease of organics with decreasing solar distance (probably under evaporation processes) and that a significant fraction of the interplanetary dust is of cometary origin [START_REF] Lasue | Inferring the interplanetary dust properties from remote observations and simulations[END_REF]. From numerous observations of 1P/Halley and C/1995 O1 Hale-Bopp in various colours, size distributions can be derived, together with the relative bulk contributions of silicate and more absorbing organic material, which remains between 33% and 60% in mass [START_REF] Lasue | Cometary dust properties retrieved from polarization observations, Application to C/1995 O1 Hale-Bopp and 1P/Halley[END_REF]. Changes noticed from innermost comae to outer comae could be attributed EPSC Abstracts Vol. 9, EPSC2014-866- to an evaporation of ices or alteration of organics, after dust is released by a nucleus.

Perspectives

Improved links between remote polarimetric observations of cometary dust and in situ data are expected from the ESA Rosetta rendezvous mission to comet 67P/Churyumov-Gerasimenko (mostly from dust instruments, i.e., GIADA, MIDAS, COSIMA) and simultaneous remote polarimetric observations (from HST and other observatories). An improved resolution in maps of interplanetary dust intensity and polarization is expected from Eye-Sat nanosatellite. This cubesat, mainly developed by students working at CNES, is expected to be launched in piggy-back on a Sun-Synchronous Orbit in 2016. Such approaches should provide better insights on solar system early evolution, including dust accretion during its formation and delivery of carbonaceous compounds on telluric planets though interplanetary dust particles impacts during its early evolution.

Acknowledgements

Support from CNES, together with the contribution of students working on the Eye-Sat project, is warmly acknowledged.


Erwan Drezen

Guillaume Rizk

Rayan Chikhi

Charles Deltel

Claire Lemaitre

Pierre Peterlongo

Dominique Lavenier

GATB: a software toolbox for genome assembly and analysis

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

DiscoSNP

Discover Single Nucleotide Polymorphism (SNP) from nonassembled reads

TakeABreak

Detects inversion breakpoints without a reference genome by looking for fixed size topological patterns in the de Bruijn graph

GATB helps you as a NGS user

The GATB C++ library gives you the opportunity to quickly develop new NGS tools that fit your needs.

GATB helps you as a NGS developer

Here is a typical workflow when working with GATB GATB-CORE transforms the reads into a de Bruijn graph, saves it in a HDF5 file that can be opened by other tools developed with the GATB-CORE API.

The core data structure of GATB is a de Bruijn graph that encodes the main information from the sequencing reads.

Strength of GATB

GATB makes this graph compact by using a Bloom filter (a space efficient probabilistic data structure) and by using a CFP additional structure that avoids false positive answers from the Bloom filter due to its probabilistic nature. 2. GATB-TOOLS: a set of elementary NGS tools mainly built upon the GATB library (k-mer counter, contiger, scaffolder, variant detection, etc.).

3. GATB-PIPELINE: a set of NGS pipeline that links together tools from the previous layer.

License & Web Site

GATB is released under the GNU Affero General Public License.

Proprietary licencing for software editors or services providers is currently being studied.

For more details on GATB: http://gatb.inria.fr

►

  de Bruijn graph: a basis for families of tools ► Data error correction ► Assembly Biological motif detection Several tools based on GATB are already available Bloocoo K-mer spectrum based read error corrector for large datasets Minia Short read assembler based on a de Bruijn graph. Results are of similar contiguity and accuracy to other de Bruijn assemblers (e.g. Velvet)

  Major facts about the GATB C++ library ► Object oriented design ► Simple and powerful de Bruijn graph API ► Simple and powerful multithreading model ► HDF5 usage for data storage ► Full Doxygen documentation including many code samples ► Test suite that checks the correct behavior of the library

  The GATB philosophy proposes a 3-layer construction to analyze NGS datasets terabytes of data. Efficient and fast NGS algorithms are essential to analyze them.ObjectiveThe Genome Assembly Tool Box (GATB)► is an open-source software ► provides an easy way to develop efficient and fast NGS tools ► is based on data structure with a very low memory footprint ► allows complex genomes to be processed on desktop computers Erwan Drézen 1 , Guillaume Rizk 1 , Rayan Chikhi 2 , Charles Deltel 1 , Claire Lemaitre 1 , Pierre Peterlongo 1 and Dominique Lavenier 1 1 INRIA/IRISA/GenScale, Campus de Beaulieu, 35042 Rennes cedex 2 Department of Computer Science and Engineering, Pennsylvania State University, USA -CORE: a C++ library holding all the services needed for developing software dedicated to NGS data.

  , D. Lavenier, R. Chikhi, DSK: k-mer counting with very low memory usage, Bioinformatics, 2013 Mar 1;29(5):652-3 R. Chikhi, G. Rizk. Space-efficient and exact de Bruijn graph representation based on a Bloom filter, Algorithms for Molecular Biology 2013, 8:22 G. Collet, G. Rizk, R. Chikhi, D. Lavenier, Minia on Raspberry Pi, assembling a 100 Mbp genome on a Credit Card Sized Computer, Poster at the JOBIM conference, 2013 Jul 1-4 (Toulouse) Best poster award. K.l Salikhov, G. Sacomoto, G. Kucherov, Using Cascading Bloom Filters to Improve the Memory Usage for de Brujin Graphs, Algorithms in Bioinformatics, Lecture Notes in Computer Science, Volume 8126, 2013, pp 364-376

a whole human genome sequencing reads can be handled with 5 GBytes of memory

	>read 1 ACGACGACGTAGACGACTAGCA
	AAACTACGATCGACTAT
	>read 2 ACTACTACGATCGATGGTCGCG
	CTGCTCGCTCTCTCGCT
	...
	>read 100.000.000 TCTCCTAGCGCGGCGTATACGC
	TCGCTAGCTACGTAGCT
	...
	How to How to analyze complex analyze complex genomes on a genomes on a simple desktop simple desktop computer ? computer ?


GATB's de Bruijn graph: a basis for families of tools ► Data error correction ► Assembly ► Biological motif detection Several tools based on GATB are already available Bloocoo K-mer spectrum based read error corrector for large datasets Minia Short read assembler based on a de Bruijn graph. Results are of similar contiguity and accuracy to other de Bruijn assemblers (e.g. Velvet)

DiscoSNP

Discover Single Nucleotide Polymorphism (SNP) from nonassembled reads

TakeABreak

Detects inversion breakpoints without a reference genome by looking for fixed size topological patterns in the de Bruijn graph

GATB helps you as a NGS user

The GATB C++ library gives you the opportunity to quickly develop new NGS tools that fit your needs.

GATB helps you as a NGS developer

Here is a typical workflow when working with GATB GATB-CORE transforms the reads into a de Bruijn graph, saves it in a HDF5 file that can be opened by other tools developed with the GATB-CORE API.

The core data structure of GATB is a de Bruijn graph that encodes the main information from the sequencing reads.

Strength of GATB

GATB makes this graph compact by using a Bloom filter (a space efficient probabilistic data structure) and by using a CFP additional structure that avoids false positive answers from the Bloom filter due to its probabilistic nature. 2. GATB-TOOLS: a set of elementary NGS tools mainly built upon the GATB library (k-mer counter, contiger, scaffolder, variant detection, etc.).

3. GATB-PIPELINE: a set of NGS pipeline that links together tools from the previous layer.

License & Web Site

GATB is released under the GNU Affero General Public License.

Proprietary licencing for software editors or services providers is currently being studied.

For more details on GATB:

  and Analysis in the Future Meeting, Santa Fé, USA, May 2014 Sequencing, Finishing and Analysis in the Future Meeting, Santa Fé, USA, May 2014 1. What is GATB ? Motivation NGS technologies produce terabytes of data. Efficient and fast NGS algorithms are essential to analyze them.ObjectiveThe Genome Assembly Tool Box (GATB)► is an open-source software ► provides an easy way to develop efficient and fast NGS tools ► is based on data structure with a very low memory footprint ► allows complex genomes to be processed on desktop computers Erwan Drézen 1 , Guillaume Rizk 1 , Rayan Chikhi 2 , Charles Deltel 1 , Claire Lemaitre 1 , Pierre Peterlongo 1 and Dominique Lavenier 1 1 INRIA/IRISA/GenScale, Campus de Beaulieu, 35042 Rennes cedex 2 Department of Computer Science and Engineering, Pennsylvania State University, USA -CORE: a C++ library holding all the services needed for developing software dedicated to NGS data.

  , D. Lavenier, R. Chikhi, DSK: k-mer counting with very low memory usage, Bioinformatics, 2013 Mar 1;29(5):652-3 R. Chikhi, G. Rizk. Space-efcient and exact de Bruijn graph representation based on a Bloom filter, Algorithms for Molecular Biology 2013, 8:22 G. Collet, G. Rizk, R. Chikhi, D. Lavenier, Minia on Raspberry Pi, assembling a 100 Mbp genome on a Credit Card Sized Computer, Poster at the JOBIM conference, 2013 Jul 1-4 (Toulouse) Best poster award. K.l Salikhov, G. Sacomoto, G. Kucherov, Using Cascading Bloom Filters to Improve the Memory Usage for de Brujin Graphs, Algorithms in Bioinformatics, Lecture Notes in Computer Science, Volume 8126, 2013, pp 364-376

http://gatb.inria.fr a whole human genome sequencing reads can be handled with 5 GBytes of memory

	>read 1
	ACGACGACGTAGACGACTAGCA
	AAACTACGATCGACTAT
	>read 2
	ACTACTACGATCGATGGTCGCG
	CTGCTCGCTCTCTCGCT
	...
	>read 100.000.000
	TCTCCTAGCGCGGCGTATACGC
	TCGCTAGCTACGTAGCT
	...
	How to How to Analyze Complex Analyze Complex Genomes on a Genomes on a Simple Desktop Simple Desktop Computer ? Computer ?


Axel Legay

Sean Sedwards

Louis-Marie Traonouez

Lightweight Verification of Markov Decision Processes with Rewards

Markov decision processes are useful models of concurrency optimisation problems, but are often intractable for exhaustive verification methods. Recent work has introduced lightweight approximative techniques that sample directly from scheduler space, bringing the prospect of scalable alternatives to standard numerical algorithms. The focus so far has been on optimising the probability of a property, but many problems require quantitative analysis of rewards. In this work we therefore present lightweight verification algorithms to optimise the rewards of Markov decision processes. We provide the statistical confidence bounds that this necessitates and demonstrate our approach on standard case studies.

Introduction

Markov decision processes (MDP) describe systems that interleave nondeterministic actions and probabilistic transitions. Such systems may be seen as comprising probabilistic subsystems whose transitions depend on the states of the other subsystems, while the order in which concurrently enabled transitions execute is nondeterministic. This order may radically affect the system's behaviour and can be optimised. By assigning numerical rewards or costs to execution traces, MDPs have proven useful in many real optimisation problems [START_REF] White | A survey of applications of Markov decision processes[END_REF]. More recently, in the context of formal verification, logics have been extended to allow model checkers to consider rewards [START_REF] Kwiatkowska | Stochastic model checking[END_REF].

Fig. 1 shows a typical fragment of an MDP. Its execution semantics are as follows. In a given state (s 0 ), an action (a 1 , a 2 , . . . ) is chosen nondeterministically to select a distribution of probabilistic transitions (p 1 , p 2 , . . . or p 3 , p 4 , etc). A probabilistic choice is then made to select the next state (s 1 , s 2 , s 3 , s 4 , . . . ). In the classic case, rewards are assigned to actions [START_REF] Bellman | Dynamic Programming[END_REF][START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF][START_REF] Baier | Principles of model checking[END_REF]. In the context of model checking, rewards are often assigned to states or transitions between states [START_REF] Kwiatkowska | Stochastic model checking[END_REF]. In both cases the rewards are summed over the length of a trace and averaged over all traces, giving the expected reward. The mechanism of accumulating rewards is unimportant to our algorithms and we need only assume that a total reward is assigned to a finite trace.

In this work we focus on MDPs in the context of statistical model checking (SMC). Model checking is an automatic technique to verify that a system satisfies a property specified in temporal logic [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]. Probabilistic (numerical) model checking quantifies the probability that a probabilistic system will satisfy a property [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF]. Statistical model checking describes a collection of Monte Carlo methods that approximate the results of probabilistic model checking.

Numerical model checking algorithms to solve purely probabilistic systems scale exponentially with the number of interacting variables in the model [START_REF] Bianco | Model checking of probabilistic and nondeterministic systems[END_REF]. Numerical algorithms to find extremal schedulers in MDPs incur the additional cost of (effectively) considering all the possible ways that the nondeterminism might be resolved [START_REF] Kwiatkowska | Stochastic model checking[END_REF]. While memory-efficient ("lightweight") Monte Carlo techniques have been developed to address the probabilistic problem (i.e., SMC), until recently [START_REF] Legay | Scalable verification of Markov decision processes[END_REF] it has not been possible to address the nondeterministic problem in this way. Building on the techniques introduced in [START_REF] Legay | Scalable verification of Markov decision processes[END_REF], in this work we present algorithms to find schedulers that (approximately) maximise or minimise the expected reward of finite traces in Markov decision processes. We use the notion of "smart sampling" to increase the efficiency of the simple sampling strategies given in [START_REF] Legay | Scalable verification of Markov decision processes[END_REF] and present statistical confidence bounds to find the maximum or minimum reward amongst multiple estimates. To demonstrate our approach, we have implemented the algorithms in our statistical model checking platform, Plasma [15,5]1 , and apply them to a number of case studies from the literature.

Schedulers and the State Explosion Problem

Given an MDP with set of actions A, having a set of states S that induces a set of sequences of states Ω = S + , a history-dependent (general) scheduler is a function S : Ω → A. A memoryless scheduler is a function M : S → A. Intuitively, at each state in the course of an execution, a general scheduler (S) chooses an action based on the sequence of previous states and a memoryless scheduler (M) chooses an action based only on the current state. History-dependent schedulers therefore include memoryless schedulers. We assume that rewards are defined by a function R : S × Ω → Q that maps (scheduler, trace) pairs to total reward. In practice, total rewards are usually defined by the accumulation of rewards assigned to actions, states or transitions.

Fig. 2 illustrates a simple MDP for which memoryless and history-dependent schedulers can give different minimum rewards for logical property X(ψ ∧XG t φ) and t > 0. The property makes use of the linear temporal operators next (X) and globally (G). Intuitively, the property states that on the next step ψ will be true and, on the step after that, φ will be remain true for t + 1 time steps. In this example, rewards (r 0 , r 1 , r 2 ) are assigned to actions (a 0 , a 1 , a 2 , respectively). In the initial state (s 0 ), both actions a 1 and a 2 can lead to traces that satisfy the property, but subsequent actions taken in state s 0 must be a 1 . If r 1 > r 2 > 0, the minimum reward will be achieved by taking action a 2 in the initial state and a 1 whenever the execution visits s 0 thereafter. To satisfy the property, a memoryless scheduler would be forced to always take action a 1 in state s 0 and would therefore not achieve the minimum possible reward.

The principal challenge in finding optimal schedulers using numerical techniques is what has been described as the 'curse of dimensionality' [START_REF] Bellman | Dynamic Programming[END_REF] and the 'state explosion problem' [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]: the number of states of a system increases exponentially with respect to the number of interacting components and state variables. This phenomenon has led to the design of sampling algorithms that find 'near optimal' schedulers to optimise rewards in discounted MDPs [START_REF] Kearns | A sparse sampling algorithm for nearoptimal planning in large Markov decision processes[END_REF] and motivates the present work.

The state explosion problem of model checking purely probabilistic systems has been well addressed by SMC [START_REF] Younes | Probabilistic verification of discrete event systems using acceptance sampling[END_REF]. SMC uses an executable model to approximate the probability that a system satisfies a specified property by the proportion of simulation traces that individually satisfy it. The state space of the system is not constructed explicitly -states are generated on the fly during simulation -hence SMC is "lightweight" and applicable to large, possibly infinite state, systems. As a consequence, SMC and may be easily and efficiently divided on parallel computing architectures.

SMC cannot be applied to MDPs without first resolving the nondeterminism. Since nondeterministic and probabilistic choices are interleaved in an MDP, memoryless schedulers are typically of the same order of complexity as the system as a whole and may be infinite. History dependent schedulers are exponentially bigger. With the exception of [START_REF] Legay | Scalable verification of Markov decision processes[END_REF], previous attempts to apply SMC to MDPs have been "heavyweight", in that they store large data structures explicitly related to the states of the system. No previous approach has attempted to approximate the rewards based algorithms of numerical model checking.

Our Approach

Using a hash function, we implement general schedulers by mapping the concatenation of an integer seed and a trace to a pseudo-randomly chosen action. By randomly sampling seeds, we are thus able to randomly sample schedulers. To find optimal schedulers we use "smart sampling" to make efficient use of a finite simulation budget. Part of the budget is used to generate a candidate set of schedulers that has high probability of containing a near optimal scheduler. The remaining budget is used to test and refine the candidate set, according to the specific confidence bounds we define for multiple estimates. Sub-optimal schedulers are removed from the candidate set and their budget is re-allocated to good ones. Hence, the performance of good schedulers is known with increasing confidence as the algorithm proceeds, until the best is known with sufficient confidence.

Related Work

There is considerable work on sampling algorithms to optimise rewards in discounted MDPs (see, e.g., [START_REF] Chang | A survey of some simulationbased algorithms for Markov decision processes[END_REF] for a survey), however this notion of rewards is significantly different to the one commonly used in model checking [START_REF] Kwiatkowska | Stochastic model checking[END_REF]. We nevertheless describe the Kearns algorithm [START_REF] Kearns | A sparse sampling algorithm for nearoptimal planning in large Markov decision processes[END_REF], since it is the basis of a recent application of SMC to MDPs [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF].

The Kearns algorithm is the classic 'sparse sampling algorithm' for large, infinite horizon, discounted MDPs. It constructs a 'near optimal' scheduler piecewise, by approximating the best action from a current state, using a stochastic depth-first search. The algorithm can work with large, potentially infinite state MDPs because it explores a probabilistically bounded search space. This, however, is exponential in the discount. To find the action with the greatest expected reward in the current state, the algorithm recursively estimates the rewards of successive states, up to some maximum depth defined by the discount and desired error. Actions are enumerated while probabilistic choices are explored by sampling, with the number of samples set as a parameter. By iterating local exploration with probabilistic sampling, the discount guarantees that the algorithm eventually converges. The stopping criterion is when successive estimates differ by less than some error threshold.

There have been several recent attempts to apply SMC to MDPs [START_REF] Bogdoll | Partial order methods for statistical model checking and simulation[END_REF][START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF][START_REF] Henriques | Statistical model checking for Markov decision processes[END_REF][START_REF] Hartmanns | On-the-fly confluence detection for statistical model checking[END_REF][START_REF] Legay | Scalable verification of Markov decision processes[END_REF]. None address the specific problem of model checking with rewards, although [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF] makes use of the Kearns algorithm.

In [START_REF] Bogdoll | Partial order methods for statistical model checking and simulation[END_REF][START_REF] Hartmanns | On-the-fly confluence detection for statistical model checking[END_REF] the authors present on-the-fly algorithms to remove 'spurious' nondeterminism from MDPs, so that standard SMC may be used. This approach is limited to the class of models whose nondeterminism does not affect the resulting probability of a property. The algorithms therefore do not attempt to address model checking problems related to finding optimal schedulers.

In [START_REF] Henriques | Statistical model checking for Markov decision processes[END_REF] the authors present an algorithm to decide an hypothesis about the probability of an MDP property. Rewards are not considered. The algorithm generates candidate schedulers by counting state-action pairs in simulations, to iteratively improve a probabilistic scheduler. The candidates are assessed using sequential hypothesis testing and the process is repeated until an example is found or sufficient attempts have been made. If found, the example is correct, but state-actions are usually shared by many schedulers, so the improvement process may actually diverge from good schedulers.

In [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF] the authors use an adaptation of the Kearns algorithm to find a memoryless scheduler that is near optimal with respect to a discounted reward scheme. By storing information about visited states, the algorithm improves on the performance of the original Kearns algorithm, but is limited to memoryless schedulers of tractable size. The resulting scheduler induces a Markov chain whose properties may be verified with standard SMC. By considering only discounted rewards, however, [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF] does not address the standard model checking problems of MDPs with rewards.

As a plausible variant of [START_REF] Lassaigne | Approximate planning and verification for large Markov decision processes[END_REF], it may seem reasonable to use the bounded properties of SMC to restrict the length of traces in the Kearns algorithm, rather than doing this implicitly with discounted rewards. The results would then be a correct approximation of numerical model checking algorithms for rewards. Unfortunately, since the Kearns algorithm relies on enumerating actions, this approach is effectively equivalent to enumerating schedulers and is not tractable.

With the exception of the original Kearns algorithm, the above approaches use data structures whose size scales with the state space of the MDP. We therefore consider these approaches to be "heavyweight". Very recently, the elements of lightweight verification of MDPs were introduced in [START_REF] Legay | Scalable verification of Markov decision processes[END_REF]. By sampling directly from history-dependent scheduler space using only O(1) memory, [START_REF] Legay | Scalable verification of Markov decision processes[END_REF] opens up the possibility of scalable verification of MDPs on parallel lightweight computational threads. The techniques of [START_REF] Legay | Scalable verification of Markov decision processes[END_REF] form the basis of the present work.

Statistical Model Checking with Rewards

SMC algorithms work by constructing an automaton to decide whether a simulation trace ω ∈ Ω satisfies property ϕ, denoted ω |= ϕ. Since simulation traces are necessarily finite, property ϕ is typically expressed in a time bounded temporal logic (see, e.g., [START_REF] Jegourel | Importance splitting for statistical model checking rare properties[END_REF]). The automaton can then be used to estimate the probability of the property or to decide an hypothesis about the probability. In particular, the expected probability of ϕ is estimated by 1 N N i=1 1(ω i |= ϕ), where ω 1 , . . . , ω N are N simulation traces selected uniformly at random from Ω and 1(•) ∈ {0, 1} is an indicator function corresponding to the output of the automaton: it returns 1 if the trace is accepted and 0 if it is not. To bound the estimation error, it is common to use the "Chernoff" bound of [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF]. The user specifies an absolute error ε and a probability δ to define the bound P(| p -p |≥ ε) ≤ δ, where p and p are respectively the true probability and estimated probabilities. The bound is guaranteed if the number of simulations N satisfies the relation:

N ≥ (ln 2 -ln δ)/(2ε 2 ) . (1)
The same notions can be used to estimate the expected reward. Given a function R(ω) ∈ [a, b], a, b finite, that assigns a total reward to simulation trace ω, the expected reward is estimated by 1 N N i=1 R(ω i ). Since rewards may take values outside [0, 1], we use Hoeffding's generalisation of (1) [START_REF] Hoeffding | Probability inequalities for sums of bounded random variables[END_REF] to bound the errors. To guarantee P(| r -r |≥ ε) ≤ δ, where r and r are respectively the true and estimated values of expected reward, N is required to satisfy the relation:

N ≥ ln(2/δ) × (a -b) 2 /(2ε 2 ) . (2)
The values of a and b are usually unknown, hence (2) cannot be implemented explicitly. We see, however, that N depends on the ratio of the absolute error ε to the range of values (a -b). The confidence of estimates of rewards may therefore be specified as a percentage. For simplicity, our algorithms use (1), assuming ε expresses a percentage as a fraction of 1.

The rewards properties commonly used in numerical model checking are based on an extension of the logic PCTL [START_REF] Kwiatkowska | Stochastic model checking[END_REF]. This extension defines instantaneous rewards (the average reward assigned to the k th state of all traces, denoted I k ), cumulative rewards (the average total reward accumulated up to the k th state of all traces, denoted C k ) and reachability rewards (the average accumulated reward of traces that eventually satisfy property ϕ, denoted Fϕ). Instantaneous and cumulative rewards are based on finite traces and can be immediately approximated by sampling, using (1) to bound the errors. Reachability rewards are based on unbounded F (the finally or eventually operator) and require additional consideration.

By the definition of reachability rewards [START_REF] Kwiatkowska | Stochastic model checking[END_REF], traces that do not satisfy the property are assigned infinite reward. Hence, if P(Fϕ) < 1 there exists a path that does not satisfy ϕ and the average reward is infinite. If P(Fϕ) = 1, all considered paths are finite and the average reward is also finite. Using sampling, it is not possible to say with certainty whether P(Fϕ) = 1, even if every observed trace of finitely many satisfies ϕ. Hence, the random variable from which samples are drawn could include the value infinity, giving it infinite variance. Error bounds that rely on finite variance, such as (1), ( 2) and the standard confidence interval, therefore cannot be applied directly.

Our solution is to implement Fϕ as F k ϕ, i.e., bounded reachability, with bound k set much longer than it is supposed necessary to satisfy ϕ. If all observed traces satisfy ϕ, the estimated average reward is within the limits defined by (1), given the hypothesis P(Fϕ) = 1 is true. If no counterexamples to F k ϕ are observed, the confidence of this hypothesis increases with increasing numbers of simulations and may be quantified according to standard hypothesis tests. In practice, the simulations required to satisfy (1) typically lead to much greater confidence with respect to the hypothesis P(Fϕ) = 1. If a counterexample is observed, the hypothesis may nevertheless be accepted according to the specified confidence. In this case the user may either conclude that the average reward is infinite, accept the calculated average reward as a lower bound or increase k and try again.

To accommodate instantaneous and cumulative rewards bounded by k, we assume that all simulation traces reach the k th state. This may be achieved by adding self loops to halting states, suitably augmented with rewards, as is standard practice in numerical model checking.

Lightweight Verification of MDPs

To avoid storing schedulers as explicit mappings, we construct schedulers on the fly using uniform pseudo-random number generators (PRNG) that are initialised by a seed and iterated to generate the next pseudo-random value. In general, such PRNGs aim to ensure that arbitrary subsets of sequences of iterates are uniformly distributed and that consecutive iterates are statistically independent. PRNGs are commonly used to implement the uniform probabilistic scheduler, which chooses actions uniformly at random and thus explores all possible combinations of nondeterministic choices. Executing such an implementation twice with the same seed will produce identical traces. Executing the implementation with a different seed will produce an unrelated set of nondeterministic and probabilistic choices. It is therefore not possible to estimate the probability of a property under a specific scheduler from either M or S.

Our solution is to use independent PRNGs to resolve nondeterministic and probabilistic choices. It is then possible to generate multiple probabilistic simulation traces per scheduler by fixing the seed of the PRNG for nondeterministic choices, while choosing random seeds for a separate PRNG for probabilistic choices. In a naive implementation, the sequence of iterates from the PRNG used for nondeterministic choices will be the same for all instantiations of the PRNG used for probabilistic choices, hence the i th iterate of the PRNG for nondeterministic choices will always be the same, regardless of the state arrived at by the previous probabilistic choices. The i th chosen action will be neither state nor trace dependent. To span the full range of general schedulers we therefore construct a per-step PRNG seed that is a hash of the integer identifying a specific scheduler concatenated with an integer representing the sequence of states up to the present.

General Schedulers Using Hash Functions

We assume that a state of an MDP is an assignment of values to a vector of system variables v i , i ∈ {1, . . . , n}. Each v i is represented by a number of bits b i , typically corresponding to a primitive data type (int, float, double, etc). The state can thus be represented by the concatenation of the bits of the system variables, such that a sequence of states may be represented by the concatenation of the bits of all the states. Without loss of generality, we interpret such a sequence of states as an integer of n i=1 b i bits, denoted s, and refer to this in general as the trace vector. A scheduler is denoted by an integer σ, which is concatenated to s (denoted σ : s) to uniquely identify a trace and a scheduler. Our approach is to generate a hash code h = H(σ : s) and to use h as the seed of a PRNG that resolves the next nondeterministic choice.

The hash function H thus maps σ : s to a seed that is deterministically dependent on the trace and the scheduler. The PRNG maps the seed to a value that is uniformly distributed but nevertheless deterministically dependent on the trace and the scheduler. In this way we can approximate the scheduler functions S and M described in Section 1.1. Importantly, the technique only relies on the standard properties of hash functions and PRNGs. Algorithm 1 is the basic simulation function used by our algorithms.

Algorithm 1: Simulate

Input:

M: an MDP with initial state s0 ϕ: a bounded temporal logic property σ: an integer identifying a scheduler Output:

ω: a simulation trace

An Efficient Iterative Hash Function

To implement our approach, we use an efficient hash function that constructs seeds incrementally. The function is based on modular division [START_REF] Knuth | The Art of Computer Programming[END_REF]Ch. 6], such that h = (σ : s) mod m, where m is a suitably large prime. Since s is a concatenation of states, it is usually very much larger than the maximum size of integers supported as primitive data types. Hence, to generate h we use Horner's method [START_REF] Horner | A new method of solving numerical equations of all orders, by continuous approximation[END_REF][18, Ch. 4]: we set h 0 = σ and find h ≡ h n (n as in Section 4.1) by iterating the recurrence relation

h i = (h i-1 2 bi + v i ) mod m. (3)
The size of m defines the maximum number of different hash codes. The precise value of m controls how the hash codes are distributed. To avoid collisions, a simple heuristic is that m should be a large prime not close to a power of 2 [START_REF] Cormen | Introduction to Algorithms[END_REF]Ch. 11]. Practically, it is an advantage to perform calculations using primitive data types that are native to the computational platform, so the sum in (3) should always be less than or equal to the maximum permissible value. To achieve this, given x, y, m ∈ N, we note the following congruences:

(x + y) mod m ≡ (x mod m + y mod m) mod m (4) (xy) mod m ≡ ((x mod m)(y mod m)) mod m (5)
The addition in (3) can thus be re-written in the form of (4), such that each term has a maximum value of m -1:

h i = ((h i-1 2 bi ) mod m + (v i ) mod m) mod m (6)
To prevent overflow, m must be no greater than half the maximum possible integer. Re-writing the first term of ( 6) in the form of ( 5), we see that before taking the modulus it will have a maximum value of (m -1) 2 , which will exceed the maximum possible integer. To avoid this, we take advantage of the fact that h i-1 is multiplied by a power of 2 and that m has been chosen to prevent overflow with addition. We thus apply the following recurrence relation:

(h i-1 2 j ) mod m = (h i-1 2 j-1 ) mod m + (h i-1 2 j-1 ) mod m (7)
Equation ( 7) allows our hash function to be implemented using efficient native arithmetic. Moreover, we infer from (3) that to find the hash code corresponding to the current state in a trace, we need only know the current state and the hash code from the previous step. When considering memoryless schedulers we need only know the current state.

Estimating Multiple Schedulers

Our algorithms require that we find the optimum estimates of a number of schedulers. We sample M schedulers and generate M corresponding estimates {r 1 , . . . , rM }, taking the maximum (r max ) or minimum (r min ), as required. To overcome the cumulative probability of error with the standard Chernoff bound (1) [START_REF] Legay | Scalable verification of Markov decision processes[END_REF], we specify that all estimates ri must be within ε of their respective true values r i , ensuring that any rmin , rmax ∈ {r 1 , . . . , rM } are within ε of their true value. Given that estimates ri are statistically independent, the probability that all estimates are less than their upper bound is expressed by P( M i=1 ri -

r i ≤ ε) ≥ (1 -e -2N ε 2 ) M . Hence, P( M i=1 ri -r i ≥ ε) ≤ 1 -(1 -e -2N ε 2 ) M , giving N = -ln 1 -M √ 1 -δ /(2ε 2
) for parameters M , ε and δ. This ensures that P(r min -rmin ≥ ε) ≤ δ and P(r max -r max ≥ ε) ≤ δ. To ensure the more usual stronger conditions that P(| r max -rmax |≥ ε) ≤ δ and P(| r min -rmin |≥ ε) ≤ δ, we have

N = ln 2 -ln 1 - M √ 1 -δ /(2ε 2 ) . (8)
Note that when M = 1, (8) degenerates to [START_REF] Baier | Principles of model checking[END_REF]. N scales logarithmically with M , making it tractable to consider many schedulers. The use of smart sampling, however, makes it unnecessary to test all candidate schedulers with maximum confidence.

Smart Sampling

The simple sampling strategies used in [START_REF] Legay | Scalable verification of Markov decision processes[END_REF] have the disadvantage that they allocate equal simulation budget to all schedulers, regardless of their merit. The idea of smart sampling is to maximise the probability of finding an optimal scheduler with a finite simulation budget and not waste budget estimating schedulers that are not optimal.

In general, the problem of finding optimal schedulers using sampling has two independent components: the rarity of near optimal schedulers (denoted p g ) and the rarity of the property under near optimal schedulers (denoted p g ). A near optimal scheduler is one whose reward or probability (depending on the context) is within some ε of the optimal value. If we select M schedulers uniformly at random and verify each with N simulations, the expected number of traces that satisfy the property using a near optimal scheduler is thus M p g N p g . The probability of seeing a trace that satisfies the property using a near optimal scheduler is the cumulative probability

(1 -(1 -p g ) M )(1 -(1 -p g ) N ). (9)
To maximise the chance of seeing a good scheduler with a simulation budget of N max = N M , N and M should be chosen to maximise [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF]. Then, following a sampling experiment using these values, any scheduler that produces at least one trace that satisfies ϕ becomes a candidate for further investigation. Since the values of of p g and p g are unknown a priori, it is necessary to perform an initial uninformed sampling experiment to estimate them, setting N = M = ⌈ √ N max ⌉. The results can be used to numerically optimise [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF], however an effective heuristic is to set N = ⌈1/p g ⌉, where pg is the maximum observed estimate (or minimum non-zero estimate in the case of finding minimising schedulers).

The best scheduler is found by iteratively refining the candidate set. At each step, the per-step simulation budget (N max ) is divided between the remaining candidates, simulations are performed and the average reward for each scheduler is estimated. Schedulers whose estimates fall into the "worst" quantile (lower or upper half, depending on context) are discarded. Refinement continues until estimates are known with specified confidence, according to [START_REF] Cormen | Introduction to Algorithms[END_REF]. With a periteration budget satisfying (1), the algorithm is guaranteed to terminate with a valid estimate.

Finding schedulers that optimise the rewards defined in [START_REF] Kwiatkowska | Stochastic model checking[END_REF] is simpler than finding schedulers that optimise the probability of a property. This is because the effective probability of rewards properties is always one. In the case of instantaneous and cumulative rewards, traces are not filtered with respect to a property, so all traces are accepted. In the case of reachability rewards, either all traces satisfy the property or the reward is assumed to be infinite. Hence, the case of probabilities less than one does not have to be quantified, just detected. The consequence of this, according to [START_REF] Hansson | A logic for reasoning about time and reliability[END_REF], is that the simulation budget to generate the initial candidate set can be allocated entirely to schedulers, i.e., N = 1 and M = N max .

Smart Reward Estimation Algorithm

ω σ j ← Simulate(M, ϕ, σ) 4 Q ← {(σ, q) | σ ∈ S ∧ Q ∋ q = N j=1 Rρ(σ, ω σ j )/N } 5 ∀σ ∈ S : trues (σ) ← 0 6 samples ← 0, conf ← 1, i ← 0 7 while conf > δ ∧ S = ∅ do 8 i ← i + 1 9 Mi ← |S|, Ni ← 0 10 while conf > δ ∧ Ni < ⌈Nmax/Mi⌉ do 11 Ni ← Ni + 1 12 conf ← 1 -(1 -e -2ǫ 2 N i ) M i 13 ∀σ ∈ S : ω σ N i ← Simulate(M, ϕ, σ) 14 samples ← samples + 1 15 Q ← {(σ, q) | σ ∈ S ∧ Q ∋ q = N i j=1 Rρ(σ, ω σ j )/Ni} 16 σmax ← arg max σ∈S Q(σ) 17 rmax ← Q(σmax) 18 ∀σ ∈ S, j ∈ {1, . . . , Ni} : trues (σ) = trues (σ) + 1(ω σ j |= ϕ) 19 Q ′ : {1, . . . , |S|} → S is an injective function s.t. ∀(n, σ), (n ′ , σ ′ ) ∈ Q ′ : n > n ′ =⇒ Q(σ) ≥ Q(σ ′ ) 20 S ← {σ ∈ S | σ = Q ′ (n) ∧ n ∈ {⌊|S|/2⌋, . . . , |S|}} 21 Z ← (trues(σmax) -samples p0)/ samples p0 (1 -p0) 22 if Z ≤ z(α) then 23 H0 is rejected
The reward property ρ may be of type instantaneous, cumulative or reachability, which are denoted I k ϕ, C k ϕ and F k ϕ, respectively, to unify the operation of our algorithms. In the case of I k ϕ and C k ϕ, k is user specified and ϕ is implicitly G k true. In the case of F k ϕ, ϕ is user specified and k is set as large as feasible to guarantee P(F k ϕ) = 1. The reward function R ρ : N × Ω → Q maps the seed of a scheduler and a trace to a reward, given reward property ρ.

Typically, the per-iteration budget will be such that the required confidence is reached according to (8) before the candidate set is reduced to a single element. Lines 10 to 14 allow the algorithm to quit as soon as the minimum number of simulations is reached.

To avoid unnecessary complication, the mechanisms by which the test of hypothesis P(F k ϕ) ≥ p 0 would achieve confidence level α for arbitrary probability p 0 are not given. Instead, to illustrate the basic idea, in lines 21 to 23 we include a simple normal approximation hypothesis test, based on the number of samples used to estimate the rewards. This is typically more than adequate.

Case Studies

We implemented Algorithms 1 and 2 in our statistical model checking platform, Plasma 1 , and were thus able to take advantage of its distributed verification algorithm on the Igrida parallel computational grid2 . All timings are based on 64 simulation cores. The following results demonstrate typical performance on a selection of standard case studies and reward properties taken from the numerical model checking literature. On these and many other examples not shown, we were able to achieve accurate results with a relatively modest per-iteration simulation budget of 10 5 simulations and a Chernoff bound of ε = δ = 0.01. The models can be found illustrated in detail on the Prism case studies website3 . Network Virus Infection Our network virus infection case study is based on [START_REF] Kwiatkowska | Probabilistic mobile ambients[END_REF] and initially comprises the following sets of linked nodes: a set containing one node infected by a virus, a set with no infected nodes and a set of uninfected barrier nodes which divides the first two sets. A virus chooses which node to infect nondeterministically. A node detects a virus probabilistically and we vary this probability as a parameter for barrier nodes. Fig. 3 illustrates the results of estimating the maximum and minimum expected number of detected attacks before a particular node is infected. Each point required approximately 15 seconds of simulation time. In Figs. 3 and4 the black lines show the true values, with their width indicating the specified ±1% error.

Self Stabilisation

The self-stabilising protocol of [START_REF] Israeli | Token management schemes and random walks yield selfstabilizating mutual exclusion[END_REF] works asynchronously to ensure that a number of networked processes share a single 'privileged status' token fairly. The protocol is designed to reach this dynamical state even if initially there are several tokens in the ring. For models containing 11 and 12 processes, we estimated the expected number of steps to reach stability, given different initial numbers of tokens. Fig. 4 plots the maximum values for 12 processes and the minimum values for 11 processes. Individual estimates required between 1 and 3 minutes of simulation time.

Gossip Protocol

The gossip protocol of [START_REF] Kwiatkowska | Analysis of a gossip protocol in PRISM[END_REF] uses local connectivity to propagate information globally. Our algorithms correctly estimate the expected minimum and maximum number of rounds necessary for the network to become connected to be 1.486 and 4.5. The correct values are 1.5 and 4.5. The average simulation time per estimate was approximately 1 minute.

Choice Coordination

To demonstrate the scalability of our approach, we consider instances of the choice coordination model of [START_REF] Ndukwu | An expectation transformer approach to predicate abstraction and data independence for probabilistic programs[END_REF] with BOUND = 100. This value makes most of the models intractable to numerical model checking, however it is possible to infer the correct values of rewards from tractable instances. The property gives the expected minimum number of rounds necessary for a group of tourists to meet. The following table summarises the results:

Prospects and Challenges

In this work we have focused on estimating the values of optimal rewards, but our techniques are immediately extensible to sequential hypothesis testing. One advantage of hypothesis testing is that schedulers which are found to satisfy the hypothesis are true examples with the specified confidence. In the case of estimation, the specified confidence is with respect to the estimate, not the optimality. Quantifying the optimality of estimates is an ongoing challenge.

Our case studies demonstrate that our approach is effective and can be efficient, but for a given hash function and PRNG, the presented algorithms typically sample from only a subset of possible schedulers. This can be addressed by also sampling from hash functions and PRNGs, however it is easy to construct examples where near optimal schedulers are vanishingly rare. This makes them difficult to find by sampling alone and motivates the development of lightweight learning techniques that accelerate convergence and draw from a larger set of schedulers. To avoid storing per-state information, or creating similar heavyweight data structures, in future work we propose to construct schedulers piecewise, by modifying and combining hash functions and PRNGs in response to information gained from simulation.

Figure 1 .Figure 2 .

 12 Figure 1. Fragment of a Markov decision process.

Algorithm 2

 2 finds schedulers that maximise rewards. The algorithm to minimise rewards follows intuitively: replace instances of 'max' with 'min' in lines 16, 17, 21 and the Output line, and replace line 20 withS ← {σ ∈ S | σ = Q ′ (n) ∧ n ∈ {1, . . . ,⌈|S|/2⌉}}. Algorithm 2: Reward Estimation Input: M: an MDP ρ ∈ {I k ϕ, C k ϕ, F k ϕ}: a reward property Rρ: the reward function for ρ H0, z(α): hypothesis P (F k ϕ) ≥ p0 and normal quantile of order α ε, δ: the reward estimation Chernoff bound Nmax > ln(2/δ)/(2ε 2 ): the per-iteration budget Output: rmax ≈ rmax, where rmax ≈ rmax and P(|rmax -rmax| ≥ ε) ≤ δ 1 N ← 1, M ← Nmax 2 S ← {M seeds chosen uniformly at random} 3 ∀σ ∈ S, ∀j ∈ {1, . . . , N } :

Figure 3 .

 3 Figure 3. Network virus infection.

Figure 4 .

 4 Figure 4. Self-stabilisation.

  rounds to converge 4.0 5.0 7.0 8.0 10.0 11.0 12.0 13.0 14.0 All the estimates are exactly correct, while the average time to generate each result was just 8 seconds.

1

  Let U prob , U nondet be uniform PRNGs with respective samples rpr, r nd 2 Let H be a hash function 3 Let s denote a state, initialised s ← s0 4 Let ω denote a trace, initialised ω ← s 5 Let s be the trace vector, initially empty 6 Select seed of U prob randomly

7 while ω |= ϕ is not decided do 8 s ← s : s 9 Set seed of U nondet to H(σ : s) 10 Iterate U nondet to generate r nd and use to resolve nondeterministic choice 11 Iterate U prob to generate rpr and use to resolve probabilistic choice 12 Set s to the next state 13 ω ← ω : s

projects.inria.fr/plasma-lab

igrida.gforge.inria.fr

www.prismmodelchecker.org/casestudies

Acknowledgement

This work was partially supported by the European Union Seventh Framework Programme under grant agreement no. 295261 (MEALS).