Boosting XML Filtering with a Scalable FPGA-based Architecture

The growing amount of XML encoded data exchanged over the Internet increases the importance of XML based publish-subscribe (pub-sub) and content based routing systems. The input in such systems typically consists of a stream of XML documents and a set of user subscriptions expressed as XML queries. The pub-sub system then filters the published documents and passes them to the subscribers. Pub-sub systems are characterized by very high input ratios, therefore the processing time is critical. In this paper we propose a"pure hardware"based solution, which utilizes XPath query blocks on FPGA to solve the filtering problem. By utilizing the high throughput that an FPGA provides for parallel processing, our approach achieves drastically better throughput than the existing software or mixed (hardware/software) architectures. The XPath queries (subscriptions) are translated to regular expressions which are then mapped to FPGA devices. By introducing stacks within the FPGA we are able to express and process a wide range of path queries very efficiently, on a scalable environment. Moreover, the fact that the parser and the filter processing are performed on the same FPGA chip, eliminates expensive communication costs (that a multi-core system would need) thus enabling very fast and efficient pipelining. Our experimental evaluation reveals more than one order of magnitude improvement compared to traditional pub/sub systems.

INTRODUCTION
Publish/subscribe applications (or simply pub-sub) are an important class of content-based dissemination systems where the message transmission is guided by the message content, rather than its destination IP address. System architectures may vary (centralized within a server or distributed over a network of brokers) but they all follow the same asynchronous event-based communication paradigm. The input is a stream of messages, generated outside of the system by a set of publishers. These messages are then selectively delivered to interested subscribers that have declared their interest by submitting profiles to the pub-sub system. This process is also known as message filtering. Examples of pubsub systems include notification websites (e.g. www.hotwire.com, news.google.com and www.ticketmaster.com), where a user can subscribe for specific events ("Rock concerts in Chicago") and get automatic notifications when the event occurs. Increasingly such environments are becoming XML-based, i.e., the messages are exchanged in XML while the users express their subscriptions using XML query languages like XPath [31].
Given the high volumes of messages and profiles, the filtering process becomes a critical performance requirement for pub-sub systems. The predominant solutions to this problem perform clustering of the user profiles based on their similarity in order to narrow down the search in the profile space. This is done by the use of Finite State Machines (FSM). In particular, elements of the user profiles are mapped to states in the state machine. The clustering is then performed by combining multiple profiles in a single FSM by analyzing and discovering the common profile paths. Since user profiles are typically known in advance (i.e., profiles play the role of data, while documents play the role of traditional queries) it is possible to be analyzed and clustered as needed before the filtering process starts.
When a document arrives in a pub-sub system, it is parsed by an event-driven parser like SAX [1] that reports low level parsing events such as: "start document", "start element", etc. As events are produced by the SAX parser, they are processed by the filtering system which uses them to drive transitions between the FSM states. For example, a transition is taken from the current FSM state if there is an outgoing edge labeled with the tag currently being processed. If during this process an "accept" FSM state is reached the document satisfies the corresponding profile(s) associated with that state.
Implementing the above approach on a traditional von Neumann architecture would requires multiple clock cycles per instruction. Consider for example, the "high level" task of identifying an "open" tag during parsing. This corresponds to multiple low level instructions (e.g. load and store) where the execution of every such instruction requires at least one clock cycle. This issue is known as the von Neumann bottleneck and can limit the filtering speed to few hundreds of clock cycles per processing a single XML tag.
Given the above bottleneck of von Neumann machines, an attempt to improve performance is to execute the tasks in parallel by adding more resources (i.e., many processors). While this idea will not eliminate the bottleneck (each processor still uses multiple clock cycles per operation) it will also create a large communication overhead between the processors. For example, one could pipeline the parsing with the filtering tasks by running them on different processors. However, when the parser produces an event it needs to notify (communicate) the filtering processor about this event (thus creating large interprocessor communication cost).
The way to resolve this limitation is to use a non-traditional highly parallel architecture. In this paper we present a novel filtering approach which is based on the use of Field-Programmable Gate Arrays (FPGA).
FPGAs are increasingly being made available as co-processors on high-performance computation systems. They are packaged in modules, which are dropped in CPU sockets on server motherboards with bridges to the FSB / Quickpath [27] [36] links on Intel platforms and Hypertransport [6] link on AMD platforms. High density FPGAs such as Xilinx Virtex-4LX 200 [37] [38] and Altera Stratix EP2S80F [2] have millions of logic gates, abundant high speed dual port memory, ALU blocks on the silicon fabric, and have high speed multi Gbps speed I/O ports [38] [3] . These high density FPGAs can be used to implement in hardware the computationally intensive portions of the software code. Multithreaded software components with streaming data input and output like the pub-sub applications are ideal candidates for acceleration on FPGA co-processing systems since a huge amount of data can be processed in parallel on the FPGA.
Since pub-sub XML filtering involves multiple queries processed over a single document data-stream, it is possible to utilize FPGAs for parallelizing the filtering performance. Each query can be implemented on the FPGA unit as a hardware datapath circuit and with appropriate optimizations it is possible to fit thousands of queries on a single FPGA chip. Moreover, having the parallel processing modules implemented on the same chip eliminates the need for expensive communications between them. This in turn allows for full pipelining of the parsing and filtering processes: as an event is produced by the parser it is immediately forwarded to the filtering module. This results in accelerated query processing and furthermore leads to substantial savings in a general purpose computation infrastructure by reducing the amount of power required by the CPUs.
In this paper we present a "proof of concept" for the use of FPGAs in boosting XML filtering performance. We utilize a four step approach that converts such query into hardware description, suitable for implementation on FPGA. The first step involves conversion of an XPath query to PERL compatible regular expressions (PCREs). The regular expressions are clustered by their common prefixes in order to produce more compact representation on the board and are then translated to VHDL using our "regex to VHDL" compiler [24]. Moreover, in order to support parent-child relationships, we introduce the use of stacks and modify the regular expression hardware to use them. The highly optimized VHDL code is then deployed on the FPGA board. The stream of documents is forwarded to the board where it is processed with high degree of parallelism. Our experimental evaluation reveals that this architecture achieves orders of magnitude improvement in the terms of running time compared to the state of the art software based XML filtering systems.
The paper is organized as follows: Section 2 presents related work. Section 3 provides in depth description of the proposed architecture. Section 4 presents an experimental evaluation of the FPGA approach compared to the state of the art software counterparts. Finally conclusions and open problems for further research appear in section 5.

RELATED WORK
One of the first works that addressed XML filtering is the XFilter [4]. This approach defines a Finite State Machine (FSM) for each XPath user profile. Every tag (element) in the profile becomes a state in the FSM, while the last tag becomes the accept state in that FSM. These machines are then executed concurrently for each message in the input. In particular, a 'start element' event drives the machine through its various transitions from state to state, while an 'end element' event makes a transition backward to a previous state. Finally, if an accepted state is reached, the document is reported as a match to the corresponding profile's subscriber. Later, the YFilter [11] system improved the matching performance by combining all profiles into a single Nondeterministic Finite Automaton (NFA). Common profile prefixes are combined and represented with a single set of states. This allows dramatic reduction in the number of states needed to represent the set of user profiles. It also improves the filtering performance of the system by processing common profile paths only once.
Other FSM-based approaches use different techniques for building the state machine as well as different types of machines. For example, [14] builds a single deterministic push down automaton using a lazy approach, [12] employs a lazily built Deterministic Finite Automata (DFA), [22] builds a transducer, which employs a DFA with a set of buffers, and [28] employs a hierarchical organization of push down transducers with buffers.
All these solutions are similar in the sense that they traverse the provided input document in a top-down fashion (i.e. in-order traversal) while advancing the state machine for each XML element (or attribute) read. Another proposed approach is to use a bottom-up traversal of the document. This idea takes into consideration the fact that an XML document typically has its more selective elements located at its leaves and uses them to perform early pruning in the query space. Examples of systems which utilize the bottomup approach include FiST [16] and BUFF [25].
The NFA based approaches discussed above are entirely software based solutions using the standard von Neumann organization. None of them takes advantage of specialized architectures to overcome the bottleneck problem which appears during XML document filtering.
Previous works [23,19,33] that have used FPGAs for processing XML documents have mainly dealt with the problem of XML parsing which in turn is transformed to implementing regular expressions on FPGAs. In particular, [23] proposes the ZuXA engine to parse XML documents. This engine employs state machines for efficient parsing based on set of rules. The paper however does not provide any discussion how this engine can be adapted to work with the XPath or twig profiles common in the pub-sub systems. A related FPGA based regular expression language parser adapted for content based routing of an XML stream has been demonstrated in [26].
There is also a large amount of research related to implementing regular expressions on FPGAs [32,18]. Here we build on our previous works [24] where we compiled PERL Compatible Regular Expressions (PCRE) to VHDL for accelerating intrusion-detection system rules using FPGAs. However, XPath query evaluation is more complex than plain regular expressions. To this end we introduce appropriate stacks that are implemented on the FPGA device.
The works in [33,19] propose the use of a mixed hardware/software architecture to solve simple XPath queries having only parent-child axis. A finite state machine implemented in FPGAs is facilitated to parse the XML document and to provide partial evaluation of XPath predicates. The results are then reported to the software part for further processing. Similarly to the ZuXA engine, this architecture can only support simple XPath queries with only parent-child axis.
There are also approaches that use specialized parallel architectures for XML processing [17,20,21]. In particular, [17] uses the Cell Broadband Engine multi-processor which consists of 8 independent processors (SPEs) that run the same software. This approach achieves parallelism by parsing (eight) XML documents in parallel at a time. Each processor implements the FSM of the ZuXA engine [23]. In addition to be only suitable for XML parsing, this solution is a combination of hardware-software approach. Similarly, the work in [20,21] addresses ways to load-balance parallel threads for low-level XML processing (e.g., XML parsing). There is also work on running XML queries over documents that are fragmented among many processors [8,10] and achieving parallelism through partial query evaluation; nevertheless, this is an orthogonal problem to filtering.
To the best of our knowledge our system is the first one to provide an entirely hardware solution to the XML filtering problem in pub-sub systems. It is also the first one able to efficiently evaluate complex XPath queries with different types of navigation directions (parent-child "/" as well as ancestor-descendant "//" axis) over the stream of XML documents. While parallelism can be achieved with multi-core machines (as a software-hardware solution), FPGAs offer a viable alternative due to their power efficiency (less power consumption and cooling costs) [34,15] as well as higher throughput. The work in [13], quantitatively demonstrates the benefits of using FPGAs over general purpose CPUs for streaming applications. While multi-core systems come with 2 and 4 CPUs it is not always feasible to achieve proportional speed-up due to the bottleneck in shared cache memory and the front side bus.

IMPLEMENTING XPATH PROFILES ON FPGAS
We start this section with a short description of the FPGA architecture and the properties that make it appealing for XML filtering. This is followed by a general overview ( Figure 2) of our compilation workflow, which loads the filtering logic on the FPGA chip. Finally we present a detailed description of the individual steps in the workflow; this includes two optimizations, namely the common prefix optimization and the area efficient character decoder.
A Field-Programmable Gate Array (or FPGA) is a semiconductor device containing programmable logic components termed "Configurable Logic Blocks" (CLB) connected trough programmable interconnections. An illustration of a typical FPGA architecture appears in Figure 1. The interconnections inside the device allow logic blocks to be interconnected as needed by the user in order to implement specific logic. Such devices allow implementation of multiple datapaths operating in parallel which makes them suitable for streaming applications like XML parsing and filtering. Moreover, because the datapath is implemented in hardware, the load and store operations from the von Neumann model are eliminated resulting in more efficient processing.  As it can be seen from Figure 2 in the first step of the compilation workflow the tag elements in the XPath expressions, representing the user profiles, are replaced with fixed length string encodings. This is done to simplify the processing and to ensure that each tag element occupies the minimum amount of area possible on the FPGA device. Reducing the footprint of the individual XML tags results in higher query density on the chip and thus better usage of the hardware.
After this step the XPath expressions are translated to their equivalent PCRE form. During this translation process the navigation directions inside the XPath expression ( parent-child "/" and ancestordescendant "//" ) are replaced with their PCRE counterparts. We describe this process in detail later in this section. In order to further reduce the query footprint on the FPGA device we cluster the regular expressions by their common prefixes. Those common prefixes are implemented as a single block on the FPGA unit. The result from the clustering step is a forest of "common prefix" trees. Each tree is compiled to generate a set of VHDL hardware blocks. The rest of the workflow involves FPGA specific compilation steps which will be discussed later as well.

Dictionary Replacement
The area of the FPGA chip is a limited resource. In order to get better usage, we minimize the tag footprint on the chip through a dictionary replacement process which replaces the XML tags in the input documents and the user profiles with fixed length strings. In

XPath to Stack-enhanced Regular Expressions
If the XPath expression contains only the ancestor-descendant axis the translation to regular expression is straightforward. While the YFilter approach, maps an XPath profile to a sequence of NFA states connected with transitions, our approach maps an XPath profile to a regular expression. As an example the XPath profile "a0//b0" will be translated to the regex " <a0> [\w\s] + [<\c\d> | </\c\d>] * <b0> ". The various regular expression operators are explained in Table 1.
The regular expression in the above example accepts a sequence of XML tags which starts with <a0> and includes <b0>. It first matches the tag <a0>. Once this is matched, it will look for one (or more) characters (the [\w\s] + part) corresponding to text between tags and then will check for any number (0 or more) of open OR closed tags (the [<\c\d> | </\c\d>] * part) before it matches <b0>.
Moreover, in order for <b0> to be a descendant of <a0> in the document, the regular expression should match before the closing of <a0>. To implement this, during the hardware generation step for this regular expression, our compiler automatically adds a negation block on </a0> so that <b0> is matched before </a0> appears in the stream. The block diagram of the regular expression on the FPGA is shown in Figure 3. Each block represents a tag parser that searches for the given tag in the document stream. The right most hardware block (depicted as a circle), provides the final result from the matching process of the regular expression. Each block receives input from the 8 bit streaming XML interface and works in parallel with the other blocks.
The translation of the parent-child axis to a regular expression requires special treatment. This is due to the fact that the regular expressions are memoryless structures and one needs to ensure that the matched XML tags occur on consecutive levels in the document. For example, the level on which the parent is matched should be remembered so as to ensure that the child is matched on a consecutive level (e.g. it is immediately below the parent). The regular expression hardware is thus modified to include the notion of memory. In our implementation this is accomplished through the use of a tag stack which keeps the current path in the XML document.
When an open tag is encountered it is pushed into the stack. Similarly when a close tag is reached it is popped from the top of the stack (TOS).
When testing a parent-child relationship, in addition to checking for the ancestor-descendant property we have to ensure that the level difference between the respective tags is one. Hence we use an extra hardware block -the TOS matching -, which continuously monitors the top of the stack and ascertains that the matched element <b0> is indeed a child of the previously matched element <a0>. In a similar way closing tags trigger the pop function and remove the head of the stack. A difference with the previous ancestordescendant match is that finding <b0> after <a0> is not enough; we need also that the top of the stack is <a0> (when <b0> is found). Since many regular expressions are using the same XML input data stream, we need only one stack block per data stream.

Common Prefix Optimization
The regular expressions derived from the XPath profiles typically depict large commonality in their prefixes. For example "a0//b0// c0//d0" and "a0//b0//c0//e0" share the common prefix "a0//b0//c0", with corresponding suffixes "d0" and "e0". The hardware cost of implementing the regular expressions is measured in terms of the FPGA area used to implement the logic. It is thus advantageous to combine multiple regular expressions into a common prefix tree. Such a tree can help reduce the area cost of the hardware by implementing the common prefix as a single block on the chip. In the above example, instead of implementing two regular expression hardware blocks and duplicating the "a0//b0//c0" logic, we can have a single implementation for the common path. As a result, more profiles can fit in a given FPGA area.
Given a set of XPath profiles, we first create their regular expressions and then sort them in alphabetical order. We then run a common prefix discovery algorithm on the sorted list of the regular expressions. The algorithm recursively grows the common prefix one tag at a time. The result is a forest of common prefix trees, each representing a set of profiles. From these trees we then create the FPGA hardware.

Area Efficient Character Decoder Hardware
Implementing XPath profiles on FPGAs mainly involves implementing character matching blocks to identify XML tags in the input document stream. The character matching hardware block [9] compares sequences of characters from the input stream to a given sequence that defines an XML tag. Figure 6 exemplifies the comparator hardware that matches an XML tag. Each character requires an 8-bit comparator block. The implemented character matching blocks for the XML tags consist of many redundant blocks, the prime examples being the open "<", close ">", and end tag "/" characters.
It is possible to simplify the character match hardware with a 8bit ASCII character pre-decoder. The character pre-decoding hardware decodes the incoming ASCII data stream at the input. An 8bit input is decoded into one of 256 possible 1-bit character signal every clock cycle. As an example, if the input was HEX "0x60", the output line for the character "a" would be high on that clock cycle and the rest of the other 255 outputs would be all zeros. The  character decoder hardware block simplifies character matching by replacing 8-bit character match hardware blocks with a 1-bit comparator and results in area efficient hardware. Figure 7 depicts the character pre-decoder block, and the simplified 1-bit comparator blocks for matching an XML tag. Moreover since 1-bit data lines are routed on the FPGA for each character in the XML tag, the FPGA routing overhead is reduced, which in turn leads to a design which offers faster clock speed.

FPGA Implementation of Regular Expressions
A regular expression syntax could be defined using various syntaxes such as PERL, UNIX, etc. Our implementation uses the PERL semantics. The compiler uses a modified version of the PCRE library v6.7 compilation flow. It simulates the behavior of a regular expression in VHDL, suitable for implementation on FPGA. We modified the compiler to take into account the stack directives and generate the hardware blocks to support parent-child axes.
After obtaining the VHDL sources for the user profiles, we generate additional hardware blocks including an input ASCII decoder, two output priority encoders (one each for queries with or without parent-child axes) and the tag stack. We group the VHDL sources into two sets, i.e. profiles without parent-child axes and profiles with parent-child axes. The organization of XPath expressions on the FPGA is depicted with an example in Figure 5. The four XPath profiles on the left correspond to expressions that contain parentchild axes and thus use the on-chip FPGA stack. When the streaming document matches a given profile, the output priority encoder is set to that profile.
We synthesize the generated VHDL code, using the XILINX synthesis tool to obtain the hardware netlist. The next step involves running the Place and Route tool, which report the clock frequency of the hardware design.
Our target FPGA is a Virtex-4 LX 200 device, and the target hardware is the Silicon Graphics RASC RC 100 board. In order for our FPGAs to run on this board we had to add a hardware module (RASC Core Services) which allows us to send and receive data and control the FPGA from the host system. Finally we generate the bitstreams that are loaded on the FPGA.

EXPERIMENTAL EVALUATION
This section describes our experimental setup and the obtained results when comparing the throughput of XML filtering performed on FPGAs, with respect to software based filtering solutions, i.e. the YFilter. This system is widely adopted as a software-based XML filtering approach. The software part of the experimental evaluation was executed on a Core 2 Quad 2.66 GHz with 8GB of RAM available. We choose YFilter because it uses more general approach for the XML filtering compared to other existing solutions. For example the lazy DFA presented in [12] has been shown to provide faster performance than the YFilter, but nevertheless assumes certain constraints for documents and profiles. We leave comparisons with such systems for future work.
In order to provide in depth evaluation of the performance for both the hardware and software implementations, we use the ToXGene XML document generator [7]. This tool generates XPath profile datasets for a specified DTD structure. We use the same set of profiles to test all methods described in this section.
We have generated multiple sets of profiles with varying path length, i.e. 2 Tags, 4 Tags and 6 Tags using the PathGenerator class in YFilter. The number of queries in each set varies from 16 to 1024. The streaming documents, used in the evaluation, vary in size from one to eight MBs.
During the experimental evaluation of the software approach we measure the throughput of the system (the size of the document set in megabytes provided as input divided by the time in seconds between the moment when the set enters the system to the moment those documents are filtered by the matching process).
For the hardware implementation we use the Silicon Graphics Altix 4700 [5] supercomputer system along with the RASC RC 100 [35] blade. We stream XML data stored in the memory (RAM) of the Altix system to the FPGAs placed on the RASC blade. We also stream the output of the priority encoders from the FPGA back to the Altix system. The output of the priority encoders is also continuously decoded by the host system, to filter the XPath expressions that have a match with the current document. As an output we provide the profile that is successfully matched as well as well the location of the match inside the document structure.
The speed in the hardware implementation is also measured in terms of throughput (MBytes/s). However we also measure the area occupied by the hardware design since it is considered a critical resource for FPGAs. In order to obtain a better understanding of the area/speed tradeoff which is something typical of FPGA based systems, we progressively increase the number of XPaths profiles processed on the FPGA. The total number varies from 16 up to 1024 profiles.

Area Utilization
With the first set of experiments we identify the impact of our two optimizations (i) the common prefix and the (ii) character predecoder on the area occupied on the chip. We consider four implementation scenarios: • Unoptimized Hardware (Unop): A system implementation without character decoding and with no common prefix optimization.
• Common Prefix Optimized Hardware (Com-P): A system which uses the common prefix optimization of the queries but without character pre-decoding.
• Unoptimized Hardware with Character Decoding (Unop-CharDec): A system that utilizes a character pre-decoding blocks, but without common prefix optimization.
• Common Prefix Optimized Hardware with Character Decoding (Com-P-CharDec): A system which takes advantage of both optimizations.
The results from these experiments are shown on Figure 8. The general trend which can be observed across the plots is that the occupied area increases linearly with increasing number of XPath queries for a given XPath length. The increased length of the queries have the same impact over the area. As expected the unoptimized hardware implementation is the one which consumes most area out of all implementational scenarios. Sometimes this can be prohibitively expensive. For example we were unable to implement the dataset which contains 1024 Xpath Queries with 6 tags with this approach because of space limitation on our FPGA.
In contrast the implementation which uses the common prefix optimization along with character decoder produces the most efficient area implementation of Xpath profiles. This approach is highly efficient when compared to the unoptimized hardware and in most cases provides five to seven times area improvement.

Performance Speedup
In this experimental set we compare the performance of both the hardware and software approaches. We use the same implementational scenarios discussed above with the same set of queries. The results of the comparison can be depicted in Figure 9.
In particular, there is a gradual reduction in throughput with the increase of the number of XPath profiles implemented on the FPGA. On average the unoptimized character pre-decoder based FPGA design for XPath filtering offers higher throughput than other designs. The design that almost always leads to the slowest speeds is the hardware implementation of common prefix optimized regular expressions without the use of character pre-decoder hardware.
Here we also compare with the performance of the software approach (YFilter). A common characteristic is that all FPGA based solutions are orders of magnitude (at cases 100 times) faster than YFilter. The performance of YFilter appears constant because it is limited from above by the bottleneck which appears in the traditional von Neumann architectures.

Summary
The area/speed tradeoff, typical for FPGAa, is apparent in the scenario with the implementation that uses the common prefix optimization with character pre-decoder. This approach provides the maximum area efficiency but with low throughput which for most cases is almost as low as common prefix optimized hardware. In general the common prefix optimization, even though improves area efficiency, brings down the clock speed, and thus the throughput of a FPGA based hardware filtering design. The second optimization that uses character pre-decoders offers better area/speed tradeoff.
The overall results of our experimentation lead to the conclusion that using an FPGA for parallel and efficient XPath filtering approach provides orders of magnitude throughput improvement (around 100 times for some datasets). It was observed that increasing XPath lengths decreased the speedup offered by FPGA. The same is also true about the increasing number of XPath profiles implemented on FPGA. Moreover common prefix optimized hardware, both with and without character pre-decoder provides better area utilization but lowers the system throughput. The reason is that, adding hardware complexity leads to lower clock rates on the FPGA. The unoptimized character decoder based FPGA implementation of XPath filtering offers the best area speed tradeoff.

CONCLUSIONS AND OPEN PROB-LEMS
This paper provides a preliminary implementation of XML filtering using a flexible FPGA architecture. Such filtering is limited in traditional Von Neumann architectures by the presence of a bottleneck between the CPU and memory. The execution of a single instruction requires multiple clock cycles for fetching, processing and storing the data back into memory. Using FPGAs alleviates this problem, by removing unnecessary operations and performing an instruction over the streaming data in a single clock cycle. Our experimental evaluation reveals order of magnitude (around 100 times) improvement in the performance speedup. We presented a hardware design and optimizations for efficiently handling XPath profile queries.
The idea of combining FPGAs and XML processing leads to many directions for further research. While our approach processes documents in a 'top-down' fashion, it is interesting to examine whether a 'bottom-up' solution is also possible. This means that document paths are first stacked and their leaf nodes are examined first for a match (which will be advantageous for documents whose more selective tags are at the leaves). A comparison with other software implemented XML filtering techniques (like lazy-DFA based approaches) is also interesting. One open problem is how to deal with dynamic updates (deletions and insertions) on the profile queries. A natural extension is to provide support for twig profiles. Unlike XPath profiles, which can only look for the presence of a given path inside the structure of the XML document, a twig pattern query identifies more complex structures like trees. To support twig profiles in our system we need a different approach that the architecture described in Section 3.
A straightforward solution for the twig pattern matching problem is to decompose the twig query into individual paths and process each path separately. The results from the individual paths are then joined together in a post processing step to produce the final outcome of the query. This approach (which would also work with our current XPath architecture) however requires extra processing time: first, there may be many path matches not related to the twig (false positives that need to be eliminated); second, the common sections of individual paths will be processed multiple times which is redundant.
Instead, a more promising approach is to employ holistic twig profile filtering, based on the Prüfer sequence [29] encodings of the XML document and the profiles. A Prüfer sequence was originally used in graph theory to describe a unique sequential encoding of a labeled tree. Since both the streaming XML document and the profiles represent trees, such encoding is easily attainable through tree traversals. This approach has been used in the past in software based XML filtering systems like PRIX [30] and FiST [16].
The main idea is to reduce the problem of twig matching to subsequence matching between the document and the profiles.
The reason why the Prüfer encoding of XML documents is an appealing method for identifying twig pattern matches within a document is because it captures well the document's structure. In particular, the sequence carries enough information to check parent-child and ancestor-descendant relationships within a tree structure. In particular, if a tree Q is a subgraph of another tree T then the Prüfer encoding of Q is a subsequence of the Prüfer encoding of T [30]. The reverse however is not true (i.e., we can have false positives).
Since the nature of subsequence matching leads also to sequential processing over the document (like parsing and filtering), we can then take advantage of the FPGA properties. We have currently an initial implementation of an FPGA architecture for supporting twig matching and are experimenting with approaches to efficiently eliminate false positives within the FPGA.
Another interesting future direction would be a comparison between multiple FPGAs and multi-core machines. Finally, an orthogonal problem is whether FPGAs can enable faster multi-query XML processing over an archived collection of documents.


extfieldsofstudy:
Computer Science
Biology
Mathematics

text:
Multiscale Computations on Neural Networks: From the Individual Neuron Interactions to the Macroscopic-Level Analysis

We show how the Equation-Free approach for multi-scale computations can be exploited to systematically study the dynamics of neural interactions on a random regular connected graph under a pairwise representation perspective. Using an individual-based microscopic simulator as a black box coarse-grained timestepper and with the aid of simulated annealing we compute the coarse-grained equilibrium bifurcation diagram and analyze the stability of the stationary states sidestepping the necessity of obtaining explicit closures at the macroscopic level. We also exploit the scheme to perform a rare-events analysis by estimating an effective Fokker-Planck describing the evolving probability density function of the corresponding coarse-grained observables.

Introduction
Over the past years there has been a rapid growth in biological knowledge even at the molecular/cell level concerning the physiology of neurons. Clinical studies and mathematical models have gone hand-in-hand enhancing our understanding and leading to breakthroughs in the field of neuroscience. One of the very first attempts in trying to systematically quantify the dynamics of nerve cells in mathematical terms can be traced back to 1907 when Lapicque introduced the modelling of a neuron's activity by an integrate and fire process using an electric circuit with a capacitor and resistor in parallel driven by a time-varying current (for a discussion see also Abbot, 1999). McCulloch and Pitts in 1943 introduced the concept of the artificial neuron while in 1952 Hodgkin and Huxley proposed the most celebrated model in the field describing the dynamics of action potentials. For their work in the area the two biophysicists won the Noble Prize in Physiology and Medicine in 1963. Other early benchmarks that shaped the field include the work of Rall [1959] who actually originated the development of compartmental models which incorporate spatial electro-physiological characteristics, the spatial distributed Fitzhugh-Nagumo model [Fithugh, 1961;Nagumo et al., 1962] describing the propagation of nerve pulses in the form of traveling waves, the work of Wiener and Rosenbluth [1946] and Greenberg and Hastings Cellular Automata [1978] modeling complex spatial patterns in excitable media with interconnected neurons, and models of spatially localized neural populations [Smith & Davidson, 1962;Griffith, 1963;Anninos et al., 1970;Wilson & Cowan, 1972]. A review on recent developments on the modelling of neuron dynamics can be found in [Inzikievitz, 2004;Herz et al., 2006].
The studies have proceeded to the development of detailed state-of-the-art models aspiring to approximate the complex dynamics underlying the physics and mechanisms of a wide range of problems stretching from the description of the behaviour of certain neural tissues such as the cerebral neocortex [Marr, 1970], the cortical and thalamic [Wilson & Cowan, 1973], the CA3 hippocumpal [Traub, 1983] to visual hallucinations [Ermentrout & Cowan, 1979;Bressloff et al., 2001], phase transitions in human hand movements [Haken et al., 1985] and working memory mechanisms [Durstewitz et al., 2000;Brunel & Wang, 2001;Laing & Troy, 2003;Durstewitz and Seaman, 2006] and from neurological disorders such as schizophrenia [Cohen et al., 1996;Hoffman, 1997;Rolls et al., 2008] and epilepsy dynamics [Babloyantz & Destexhe, 1986;Theoden et al., 2004] to therapeutic surgical procedures such as deep-brain stimulation process [Rubin & Terman, 2004] to name just a few.
Due to the strongly heterogeneous nature of the nerve cells as well as their stochastic and nonlinear very-large scale interactions, the collective -macroscopic -behaviour of such multicellular systems is far from intuitive [Destexhe & Contreras, 2006] and cannot, most of the time, be straightforward predicted. These largescale subunits interactions result to a remarkable and fundamental feature of analogous systems: the collectiveemergent-macroscopic -behaviour is far from intuitive [and cannot, most of the time, be straightforward predicted. Sustained oscillations, travelling waves, multiplicity of stationary states and spatio-temporal chaos [Kopel & Ermetrout, 1986;Kaneko, 1992;Chow, 1998;Freeman, 1999;Bressloff & Coombes, 1999;Chow & Koppel, 2000;Coombes & Osbaldestin, 2000;Haskell et al., 2001;Roxin et al., 2004;Coombes, 2005] are some examples of the rich nonlinear dynamics emerging in the coarse-grained population-level.
For the systematic analysis of such behaviours, one of the most critical issues in the area is related to the bridging of the different scales: the micro-scale, where the interactions of the subunits take place and the macroscopic scale where the dynamic behaviour of the constitutive system/ tissue/organ/ emerges. Traditionally, the gap between high-dimension (microscopic) and low-order (coarse-grained level) is bridged through closures, relating higher-order, fast moments of the evolving detailed distributions to a few, low-order, slow, "master" moments of the underlying detailed distributions [see for example : Treves, 1993;Bressloff et al., 1997;Bressloff, 1999;Haskell et al., 2001;Casti et al., 2002;Cai et al., 2004;Rangan & Cai, 2006]. However no general systematic methodology for deriving such closures exists; they are often based on ad hoc assumptions such as identical subunits, completely regular connections, and infinite populations, that may introduce a systematic bias and possibly miss important qualitative at the coarse-grained level [see for example the critical discussion on the influence of certain closure approximations applied to integrate and fire networks in [Ly& Tranchina, 2007].
Good closures in the form of ordinary, or partial differential equations (ODEs or PDEs), describing the coarse-grained dynamics may, in principle, be extracted from the individualistic interactions, but usually are unavailable-or extremely difficult to derive. The problem is that due to heterogeneities, nonlinearities and stochastic interactions serious limitations come up when trying to systematically bridge the different time and spatial scales, in order to extract effective-macroscopic equations.
Without the availability of such good closures, what is usually done is the exploration of the coarsegrained dynamics by performing long runs in time using the available microscopic simulator starting from different initial distributions and averaging over many ensembles. However this approach is both computational consuming and inadequate for important tasks such as the systematic detection of coarse-grained criticalities, stability analysis and the computation of coarse-grained bifurcation diagrams. There is an arsenal of state of the art computational tools that are good for such system-level analysis, yet these rely on the availability of closed form coarse-grained models.
For our illustrations we used a simple dynamic individual--based model based on the work presented in [Kozma et al., 2005]. Here, each neuron can be in one of two states: excited (activated) or non-excited (inactivated). Neurons interact with their neighbors on a random regular network [Bollobás, 2001;Newman, 2003] in a bidirectional way and change their states over discrete time in a probabilistic manner according to simple rules involving their own states and the states of their links. The relatively simple rules governing the interactions of the subunits result to the emergence of complex dynamics in the coarse-grained level including phase transitions from low to high density of excited neurons and multiplicity of coarse-grained stationary states.
The Equation-Free approach is exploited to perform system-level tasks such as bifurcation, stability analysis and estimation of mean appearance times of rare events, circumventing the need for obtaining coarsegrained models in closed form. Using the microscopic simulator as a black-box timestepper of the network excitation density we compute the coarse-grained equilibrium bifurcation diagrams and examine the stability of the solution branches with respect to a parameter describing how the state of a neuron affects the state of its links. Coupling the Equation-Free approach with Kramers' theory [Hänggi et al., 1990;Haataja et al., 2004;Sriraman et al., 2005] we also compute the mean transitions rates between the apparent coarse-grained equilibria. The required closures are obtained by relating higher-order moments of the state distribution with lower-order ones "on demand' combining the Equation-Free approach with Simulating Annealing.
The paper is organized as follows: in section 2.1 we give a short description of the basic concept of the Equation-Free framework; in section 2.2 we describe how the Equation-Free approach combined with Simulated Annealing can be exploited to enable microscopic simulations to converge on the coarse-grained slow manifold, while in section 2.3 we show how it can be combined with Kramers' theory for rare-event calculations. In section 3 we describe the individual-based model and in section 4 we present and discuss the simulation results; finally we conclude in section 5.

The basic concept
Consider a complex system for which the only good modelling description is available on the microscopic (atomistic, molecular, agent-based) level. We also believe that macroscopic laws in the form of ordinary and/or partial-integro-differential equations (ODEs/PDEs) which can describe the coarse-grained dynamics of the underlying physics may be in principle derived but due to the inherent complexity of the problem, accurate quantitative closures have not yet been extracted. For the sake of presentation, lets us write these macroscopic equations in the form of where denotes a low-dimensional coarse-grained state vector, observable through measurement, is a vector denoting the system parameters, and is a not explicitly available operator. The question is how to systematically perform system-level tasks (e.g. construct bifurcation diagrams of steady and time-dependent solutions and perform their stability analysis) when direct simulation of the microscopic code is the only possible choice and despite the fact that F is not explicitly available.
A brief schematic of the concept of the coarse timestepper is given in Figure 1. The coarse timestepper is a fundamental building block in the Equation-Free framework. It is the way to efficiently obtain macroscopic input-output information from a microscopic simulator bypassing the need to extract macroscopic equations in a closed form. Let denote the distribution function of the microscopic variables over the set of microscopic coordinates x. The main assumption behind the coarse timestepper is that a coarse-grained model for the fine-scale dynamics exists and closes in terms of a few coarse-grained variables, say, . Typically these are usually the low-order moments of the microscopically evolving distributions. The existence of a coarsegrained model implies that the higher order moments, say, , of the distribution become, relatively quickly over the time scales of interest, "slaved" to the few lower ones . This scale separation could be viewed in the form of a singularly perturbed system, reading: Eq. (2a) corresponds to the "slow" dynamics and Eq. (2b) to the "fast" ones. Notice that under certain conditions [Fenichel, 1979] the overall dynamics defined by Eq. (2) will approximate very quickly the dynamics of the following system The above system can-under some mild assumptions [Fenichel, 1979;Gear et al., 2005]-be written in the form of Eq. (1) by applying the implicit function theorem for eq. (3b), i.e., where q is a continuous function defining the relation between the slow and the fast variables after a short time horizon, i.e. the slow manifold on which the dynamics of the system evolve after a fast transient phase (see figure 2).
What a coarse time-stepper does, in fact, is providing such closures on demand (relatively short bursts of fine scale simulation naturally establish this slaving relation (refer to [Gear et al., 2002;Makeev et al., 2002;Runborg et al., 2002;Kevrekidis et al., 2003;Siettos et al., 2003;Kevrekidis et al., 2004;Gear et al., 2005] for more detailed discussions). Briefly, once the appropriate macroscopic observables have been identified, the coarse timestepper consists of the following essential components: (a) Prescribe a coarse-grained initial condition (e.g. the density of activated neurons in the network) ( ) Evolve these realizations in time using the microscopic simulator for the desired short macroscopic time T, generating the value(s) . The choice of T is associated with the (estimated) spectral gap of the linearization of the unavailable closed macroscopic equations Siettos et al., 2003].
The above procedure can be considered as an "black box" coarse timestepper that, given the initial state of the reports the solution, the result of the integration, after a given time horizon T , i.e., , where is the temporal evolution operator of the system. At this point one can: T Φ (e) Implement a computational superstructure such as Newton's method [Kelley, 1995] or Arnoldi's method [Saad, 1992;Christodoulou and Scriven, 1998] or Newton-Picard method (e.g the Recursive Projection Method of ; Kelley et al., 2004) as a shell around the coarse-timestepper to extract systems-level information. For example coarse-grained (macroscopic) steady states can be obtained as fixed points, using T as sampling time, of the mapping : If on the other hand, a periodic oscillatory behavior is observed then one seeks for solutions which satisfy with T denoting now the period of oscillation. Periodic solutions can be computed as fixed points of the mapping (5) augmented by the so-called phase constraint (also called a pinning condition) which factors out the infinite members of the family of periodic solutions [Lust et al., 1998;Kavousanakis et al., 2008].

Using the Equation-Free approach to converge on the slow manifold
For the scheme to be accurate, the overall procedure has to be applied when the system evolves on the slow manifold (i.e., the fixed point iteration (3) has to be solved on the slow manifold). However, in general, the lifting operator creates microscopic distributions off the slow manifold. If the time required for trajectories emanating from initial conditions off the slow manifold to reach the slow manifold is very small compared to T the above requirement is satisfied for any practical means. But still we can farther enhance our calculations by forcing our system to start from consistent to the coarse-grained variables (lower-order moments of the microscopic distribution) microscopic initial conditions on the slow manifold. This can be achieved by exploiting the Equation-Free approach as follows (see figure 3):
The coarse timestepper is used to approximate from short runs the drift ) (ψ u and from the time dependence of the variance the diffusion coefficient D of a Fokker-Planck equation of the form is the probability density function of the coarse-grained coordinate ψ, ) (ψ u is the drift and ) (ψ D is the diffusion coefficient. The drift and the diffusion coefficients can be estimated by appropriate initialized short-time simulations sing the coarse timestepper from t llowing relations: When the potential barrier between the stable and the unstable point is relatively big to the noise, the mean escape rates between two coarse-grained equilibria can be estimated by the formula [Erban et al., 2006] : stands for free-energy of the system, which is a function of the drift and the diffusion coefficient given by:

Simulations results and discussion
The simulation results were o ed using a total of N = 20000 neurons, while d is set equal to four. The random regular graph was generated usi the proce cribed in [Steger & Wormald, 1999]  ). As it is clearly seen, the model exhibits some interesting nonlinear dynamics. Specifically, for small values of the parameter ε and depending on the initial conditions there are two stable stationary states one corresponding to the "all-off" state (all neurons turn to be inactivated) and the other to a partially activated network (see figure 5a). As the value of ε increases -and now independently on the initial conditions -the second stable stationary point disappears and the network converges to the "all-off" state which is the only possible one (figure 5b). The "all-off" state looses its stability for bigger values of ε making the network to converge to a single non-zero stable stationary state (figure 5c, d). The above results are indicative of the existence of several critical parameter values marking the nset of qu o alitative transitions in the dynamic behavior of our system -and as the bifurcation theory suggeststhe emergence of unstable coarse-grained states which are unreachable through plain long-run temporal simulations.
At this point it is interesting to note that in the case of the two stable stationary states the non-zero is a metastable one: due to stochasticity the network's state is prone to jump to the "all-off" region of attraction (see figure 6). Based on the above simulation results, there are three obvious questions to be answered here: (a) which is the complete set of fixed point solutions (both stable and unstable ones), (b) which are the critical points which mark the onset of the phase transitions, and, (c) what is the transition rate of a stable metastable state merging f ε.
e or small values of We answer all three questions by exploiting the Equation-Free approach. The first two questions questions are being answered by performing coarse-grained bifurcation analysis while the third one by exploiting Kramers' theory for rare-events calculations.
As described in section 2.1 the first step in the Equation-Free framework is the selection of the coarse-grained statistics. The assumption is that a good closure in principle exists in terms of a set of a few coarse-grained observables which correspond to some lower order moments of the microscopic distribution. This means that the long-term dynamic behavior of the problem lies on a low-dimensional manifold parametrized by these few coarse-grained observables. Higher-order moments become quickly, relative to the macroscopic time, functions of the lower-order ones. But in what level does this happen? In order to give an answer to this question we constructed the phase portrait of the density of the activated neurons, which corresponds to the first order moment of the microscopic distribution and the densities of pairs between the states of the neurons, which correspond to the second spatial moment of the network distribution [Murrell et al., 2004].
These are the densities of the activated neurons linked to activated neurons ({11} pairs), the activated neurons linked to inactivated neurons ({10} pairs) and the inactivated neurons linked to inactivated neurons ({00} pairs). Under this pairwise representation, the density of {01} pairs is equal to the {10} one. Figure 7 depicts the phase portraits of the p vs {10} pair density around a nonzero stationary state at ε = 0.14. The trajectories were obtained by averaging over 500 ensemble realizations. We were able to create microscopic state configurations consistent with the first, second and third (which correspond to triples) -order moments by utilizing the SA procedure, thus initializing these coarse-grained quantities at will. As it is clearly shown after an initial short transient that the dynamics are evolving on a one-dimensional slow manifold, parametrized by the first-order moment-towards the coarse-grained stable stationary state; the higher order moments become unc ions o f t f the first one. This implies that the expected dynamics over the ensembles conditioned on this slow m ifol hrough SA closes and hence can be described by just the first-order moment, i.e. the density of a d t activated (or inactivated) neurons.
Hence it makes sense to search for coarse-graine bria as fixed points of the coarse timeste n d equili pper the result of lifting consistently an initial distribution conditioned on the slow manifold using SA with dT = 1 time step and evolving on the slow manifold for a time interval of T = 5 time steps to the final distribution restricted back to the first-order moment by averaging over 10000 realizations. The coarse-grained derivatives for the fixed point iteration (here, Newton-Raphson) were estimated numerically through the coarse timestepper, by appr using the coarse timestepper-based fixed point lgorithm over 10000 consistent to the coarse-grained variable p microscopic network realizations, augmented by the pseudo arc-length continuation [Chan & Keller, 1982]: opriately, perturbing the coarse-grained initial conditions, constrained again with SA on the slowmanifold.
The fixed-point solutions of (12) were obtained a ( ) The procedure involves the iterative solution of the following linearized system where k is the iteration counter and is the approximation of the solution at the kth iteration. is the jacobian matrix: which can be approximated numerically by calling the timestepper at appropriately perturbed values of the corresponding unknowns. In our difference approximation we used perturbations of the order of 10 -2 . The above computational framework enables the microscopic simulator to converge, to its own coarse-grained stable and unstable (which are not reachable through long run temporal simulations) points and their locations with respect to the bifurcation parameter. Τhe eigenvalues (here is just one) of the Jacobian d dp  As demonstrated already in figure 6, the relatively high activity stable phase is metastable: after a relatively long period of time, depending on the value of ε (the closer to the turning point the smaller and vice versa) the network can in principle "jump" from a high activation to the "all-off" phase. The mean rates for a particular value of our bifurcation parameter ε was calculated in a systematic, computational efficient way by approximating the coarse-grained free-energy of the network, as described in section three.
gure 9 depicts the free energy as computed by short runs and 100000 realizations at ε = 0.162. Here the reaction coordinate is the distance from the density of activated neurons p at the node, i.e. node p p − = ψ . Again in order to calculate the drift and the diffusion coefficient as they are defined in (9), we used the procedure of section 2.3: for each ψ we lifted to consistent -to the coarse-grained variable pmicroscopic network realizations whose second-order moments (the densities of pairs) lie on the slow manifold. Using Eq. (10) we found a mean first-passage time of the order of 20000 time steps which is in good agreement ith the outcomes of the long-run temporal simulations.
] that can be exploited to extract the right coarse-grained variables hen dealing with more complex physics. easure 8.3 of PERATIONAL PROGRAM "COMPETITIVENESS" in the 3rd Community Support Program. w

Conclusions
We have demonstrated how the Equation-Free approach, a computational framework for multi-scale computations can be exploited to systematically analyse the coarse-grained dynamics of a network of neurons. For our illustrations we used a extremely simple individual-based model which, however, is able to catch a fundamental feature of such problems which is the emergence of complex dynamics in the coarse-grained level including phase transitions from high to low (or even to zero) density of excited neurons and multiplicity of coarse-grained stationary states. We used the concept of the coarse timestepper and that of Simulated Annealing to enable the microscopic simulator to converge to its own coarse-grained slow manifold on which the coarsegrained dynamics evolve after a fast-in the macroscopic scale-time interval. Using the scheme we were able to track the coarse-grained equilibria (stable and unstable) and examine their stability with respect to a parameter representing the neurons excitation probability. We detected the coarse-grained critical points marking the onset of the phase transitions and we derived a Fokker-Planck equation to approximate the mean escape rates of the metastable coarse-grained equilibria. Further research could proceed along several directions including the systematic analysis of models closer to the real-physics, the use of more realistic network structures, advanced techniques such as the Computational Singular Perturbation method for the systematic reconstruction of the slow manifold [Lam & Goussis, 1994;Goussis & Valorani, 2006] as well as the use of Diffusion Map techniques [Coifman et al., 2005;Kolpas et al., 2007 w         the stationary state is marked with the square marker. The initial conditions corresponding to different densities of pairs on the graph were created by simulated annealing. The dynamics are characterized by two stages: a fast approach (within around 3 time steps) to a one dimensional slow manifold parametrized by the first-order moment and then a slower decay to the stationary state (marked with the square) on the "slow manifold" on which all long-term dynamics evolve. Figure 8. The one-dimensional coarse-grained bifurcation diagram of the density of activated neurons p (firstorder moment) with respect to model parameter ε. Solid lines correspond to stable coarse-grained equilibria, while, the dotted line corresponds to unstable ones. The inset depicts the computed eigenvalue determining the systems stability. Figure 9. The coarse-grained free energy for ε = 0.163 with respect to the distance, ψ from the metastable point.
The minimum corresponds to the coarse-grained metastable equilibrium and the maximum to the unstable one.


Adaptive Continuous Query Reoptimization over Data Streams

A data stream is a series of massive unbounded tuples continuously generated at a rapid rate. Continuous queries for data streams should be processed continuously, so that a strict time constraint is required. In most previous research studies, in order to guarantee this constraint, the evaluation order of join predicates in a continuous query is optimized using a greedy strategy. However, because a greedy strategy traces only the first promising plan, it often finds a suboptimal plan. To reduce the possibility of producing a suboptimal plan, in this paper, we propose an improved scheme, k-Extended Greedy Algorithm (k-EGA), that simultaneously examines a set of promising plans and reoptimize an execution plan adaptively. The number of promising plans is flexibly controlled by a user-defined range variable. The scheme verifies the performance of the current plan periodically. If the plan is no longer efficient, a newly optimized plan is generated. The performance of the proposed scheme is verified through various experiments to identify its various characteristics.

Introduction
Many recent research studies of emerging applications often need to deal with a data stream, which is a massive unbounded sequence of data elements continuously generated at a rapid rate. Such applications can be web click monitoring, sensor data processing, network traffic analysis, telephone recording and multimedia data processing. In these applications, a continuous query is generally used to monitor the ongoing information on data streams. It is registered in advance and produces results continuously and repeatedly whenever new data elements of data streams are generated. Since new data elements can be rapidly generated and should be processed in real time, a continuous query should be evaluated as quickly as possible. For this, query optimization is essential. However, query optimization techniques used in a traditional database management system are not effective for data stream applications because of the unpredictable characteristics of data streams [1]. For this reason, many previous research studies [1], [4]- [6], [8], [9] are focused on devising a method of optimizing the execution of continuous queries. A multiway join query contains multiple join operations.
Finding the optimal execution sequence of join opera- tions in such a query is problematic since a join operation is very expensive. The execution plan of a multiway join query should be continuously reoptimized in data stream environment. For this purpose, properties such as the input rates of its operand data streams and the join selectivities of its join operations should be continuously monitored. However, such monitoring requires a huge runtime overhead. In [7], given a join operation, the domain sizes of its two join attributes are monitored to estimate its join selectivity. In [5], STREAM actually obtains the join selectivity of every join operation in every possible left-deep tree plan by additionally executing them for sampled tuples. Therefore, it can obtain more accurate statistics, but it must evaluate a number of plans in run time.
To estimate the current evaluation cost of a join operation, the past properties of the operation are used. Even though the properties of a join operation may be changed in runtime dynamically, they may not vary randomly. Instead, they tend to have valid ranges respectively. This means that the past properties of a join operation can be used to predict the current cost of the operation more accurately. For this purpose, the properties of join operations that are executed at least once are maintained in a table called cost Statistics.
In order to continuously reoptimize the query, the time to generate the optimized plan of a multiway join query is also an important factor. Consequently, most approaches employ a greedy optimization strategy with respect to either the size of the result [4] or join selectivity [5]. Even though the greedy strategy can generate its solution very quickly, it is apt to produce a suboptimal plan. This is because there is no backtracking mechanism to undo the intermediate order of join operations when the order turns out to be nonoptimal. However, employing a backtracking mechanism in data stream environments is also inadequate because of the time constraint as mentioned above.
In this paper, we propose a query optimization scheme called k-Extended Greedy Algorithm (k-EGA) to produce the optimal execution plan of a multiway join query. As mentioned above, the greedy strategy traces only one promising subplan, whereas k-EGA traces a set of promising subplans simultaneously. Intuitively the cardinality of the set affects the optimality of the algorithm as well as optimization time. As it increases, the possibility of generating a more optimal plan is increased, but the time to produce the plan is increased as well. To identify the effective set of promising subplans, given a user-defined range variable k, the cost of each subplan is represented as a range based on the aver-age cost and its standard deviation in the cost statistics table. Among the join operations of a multiway join query, the one with the smallest average cost is chosen first. If the range of the selected join operation overlaps with the remaining join operations, they are traced together. As k becomes larger, the range of each join operation becomes broader. This means that the number of traced subplans increases, thereby improving the optimality of the optimized plan of a query. However, the time to produce the optimized plan clearly increases as well. By flexibly controlling the number of simultaneously traced subplans, our proposed scheme provides a desirable trade-off between the optimality of an execution plan and the time to generate the plan. Depending on the current processing load, the range variable k can be flexibly adjusted.
The remainder of this paper is organized as follows: In Sect. 2, related work is presented. In Sect. 3, the join cost model is described on the basis of the cost statistics for a multiway join query. In Sect. 4, the proposed method that can find the optimal execution plan of a multiway join query is introduced. In Sect. 5, the performance of the proposed method is analyzed through a series of experiments. Finally, In Sect. 6, our conclusions are presented.

Related Work
There are two common techniques for processing a join operator. The first technique, involving a pipelined processing tree (PPT) [17], uses a multiway nested loop to generate the result of a multiway join query without producing any intermediate join result. The second technique uses a binary linear processing tree (BLPT) [17] that gives intermediate join results. A representative research on the first technique is a multi-stream join algorithm called MJoin [1]. It produces the result of a multiway join query without any tree-structured plan. In [4], a greedy-based optimization method is proposed to determine the evaluation order of each operand of an MJoin operator.
STREAM [5] uses the A-Greedy algorithm to adaptively determine the evaluation order of an MJoin operator to minimize the processing cost of a multiway join query whose join selectivities vary unpredictably over time [14], [15]. The ongoing change in the join selectivity of every possible execution order is continuously monitored to adaptively change the evaluation order of each operand in an MJoin operator. In [18], A-Caching focuses specifically on the problem of the adaptive placement and removal of caches in an MJoin operator to optimize join performance. It provides algorithms for selecting caches, monitoring their cost and benefits in current conditions, allocating memory to caches, and adapting as conditions change.
An eddy [11] accomplishes very fine-grained adaptivity without employing any fixed execution plan. Its query execution model routes input and intermediate result tuples among query operators by making independent routing decisions for each of the tuples to reflect the dynamically changing runtime selectivity of each query condition. Based on the eddy query processing framework, a continuously adaptive, continuous query (CACQ) [8] is implemented. It provides a grouped-filter index to simultaneously evaluate multiple selection predicates and a space-efficient generalization of double-pipelined hash joins called SteMs [16]. The execution order of multiple join operations in a multiway join query is dynamically determined by a ticket-based routing policy that executes a lottery between query operators. The winning chances of a particular operator are proportional to the number of tickets it owns. The operator with a higher selectivity gets more tickets, which means the operator will receive more tuples early through the eddy [11].
In [7], as a modified architecture of the original eddy that uses SteMs, a query operator called a STAIR (Storage, Transformation and Access for Intermediate Results) is introduced to perform a join operation. A STAIR holds the runtime state of join operations. It employs the second technique, BLPT. Even though it solves the drawbacks of the original eddy, they address the problems caused by allowing the query state to be modified and migrated across STAIRs during query execution. STAIRs also solve the problem in which the ability of an eddy to adapt is constrained by the state that accumulates in the query operators. On the basis of the cumulative selectivity of every subplan of join operations, a greedy optimization strategy is used to determine the execution plan of STAIRs dynamically.

Cost Statistics
To define the cost of a binary join operation, we use a unittime-based cost model [13] that defines the evaluation cost of a binary join operation to be the processing cost of all incoming tuples in a unit time. In the execution of a binary linear processing tree, the result of one binary join operation should be materialized because it can be one of the operand for the next binary join operation. Consequently, the result size of a binary join operation affects the overall evaluation cost of a multiway join query, so that the unit-time-based cost model is modified for a multiway join query in this paper. The terms used in the modified cost model are described in Table 1. A join operation R S should be executed when a new tuple arrives from any of its two operand streams R and S . Accordingly, the evaluation cost C i (R S ) of a binary join operation R S for newly arriving tuples in unit time should be symmetrically defined as follows.
The generating cost of the result of a binary join operation, C s (R S ) is proportional to the size of the result and is defined as.
Therefore the total evaluation cost of a binary join operation on data streams R and S becomes A multiway join query can be evaluated using a number of different execution plans. All combinations of the join operations that can appear in the plans are maintained in cost statistics. A row in cost statistics is formally defined as follows.
) and σ(C(J q m )) denote the average cost of the J q m and its standard deviation. The term count(J q m ) denotes how many times the join operation J q has been executed before. If count(J q m ) = 0, the join operation J q m has not been executed before, so that its current statistics are estimated.
When a join operation is not executed, the cost of the join operation is estimated using Eq. (3), and the estimated selectivity is reflected to its statistics in the CS table. When the join operation is part of the optimized execution plan, its real selectivity is obtained while executing the plan and is reflected onto the CS table. If this is the first real selectivity of the join operation, its past estimated statistics are deleted and its statistics are replaced with the real selectivity. In addition, its count field is set to 1. The statistics of a join operation in the CS table represent how the cost of the join operation fluctuates. If the cost of a join operation fluctuated widely in the past, its standard deviation becomes large. This is a simple three-join operation between the streams R, S , T and U. Table 2 shows the cost statistics table of the query Q when the execution plan (R S ) (T U), has been executed three times and the other plan (((S T ) R) U) has been executed five times before. The values of μ and σ are merely example values.

Subplan and Set of Traced Subplans (STSs)
A complete execution plan for a multiway join query with n join operations is a binary linear processing tree that determines the execution order of n join operations. The objective of query optimization is to find the execution order with the smallest cost. The query can be evaluated using a number of different execution plans. To denote part of an execution plan, the term m-subplan is defined formally in Definition 2. The evaluation cost of each subplan is defined by the accumulated cost of all the join operations in the subplan. Since the cost of an individual join operation is maintained by the average cost and its standard deviation in the CS table, the cost of the subplan is also represented in the same manner as follows:

Definition 3. (Cost of m-subplan)
Given a multiway join query q with n join operations J q = {J 1 , . . . J n }, its child plan sp m+1 is sp m ∪ {J e }, whose cost is obtained as follows: i) The average cost of sp m+1 , μ(sp m+1 ) = μ(sp m ) + μ(J e ). ii) Suppose that sp m and J e are independent, the standard deviation of sp m+1 , denoted by σ(sp m+1 ), is σ(sp m ) 2 + σ(J e ) 2 + 2cov(sp m , J e ). cov(sp m , J e ) is a covariance, which is a measure of how much two variables, σ(sp m ) and σ(J e ), change together. Since sp m−1 and J e are independent, cov(sp m , J e ) = 0, so that The effective range of an m-subplan is defined by a range variable k in Definition 4. The minimum value of the range, i.e., μ(sp m ) − k*σ(sp m ), is denoted byĈ min (sp m ) and the maximum value of the range, i.e., μ(sp m ) + kxσ(sp m ), isĈ max (sp m ).
As mentioned above, k-EGA simultaneously traces a set of promising subplans, which is defined by a set of traced subplans (STS) as follows.
The standard deviation of the cost of a subplan represents how much the cost of the subplan fluctuates. When the value is small, the cost is almost the same in the past, which indicates the cost can accurately represent the current cost of the subplan. In contrast, when the value is large, the cost estimation based on the average cost is not accurate enough. Our scheme traces not only the subplan with the smallest average cost, but also those subplans whose effective ranges overlap with the subplan together. The subplan whose cost cannot be estimated exactly is most likely to be traced because the effective ranges of the subplan are wide. Consequently, the subplan is examined more closely. This is because one of the subplans may lead to a better plan.
Intuitively, the range variable k affects the cardinality of its STS m . The greater the value is, the wider the effective range of each subplan becomes. Consequently, the size of STS m increases as well. Furthermore, it means that more promising subplans are traced simultaneously to find out the final optimized execution plan. Therefore, the optimality of the final plan would be better than the one obtained by the general greedy method. On the other hand, the time to generate the final plan may be increased. The proposed method provides a good trade-off between the optimality of the final plan and the time to generate the final plan since they can be flexibly controlled by the range variable k. Figure 1 illustrates how a range variable k plays a role in controlling the cardinality of STS. Suppose that there are four m-subplans {sp m 1 , sp m 2 , sp m 3 , sp m 4 } and that they are the candidate m-subplans for a multiway join query with n join operations (m ≤ n). To determine the STS m of the query, the subplan with the smallest average cost is chosen first as in the greedy optimization strategy. Consequently, the subplan sp m 3 is selected first. The effective range of each of these four subplans is calculated. Suppose that the range variable k is set to 1, then the effective range of sp m 3 be-

k-Extended Greedy Algorithm
The overall steps of k-EGA are very similar to those in the greedy optimization strategy that determines the execution order of one of unordered join operations in each optimization step. Likewise, given a multiway join query q with n join operations J q = {J 1 , . . . J n }, k-EGA determines STS 1 in the first step and then generates its child subplans STS 2 . This procedure is continued until STS n is found. Among the n-subplans in STS n , that with the minimum evaluation cost sp n min (sp n min ∈ STS n ) is the final optimized plan. Given STS m = {sp m 1 , sp m 2 , .....sp m k }, STS m+1 is generated. For each m-subplan sp m i (1 ≤ i ≤ k), one additional join operation J which is not part of sp m i is added to sp m i to generate a candidate (m + 1)-subplan. Although there are (n − m) different (m + 1)-subplans for each sp m i , only those (m + 1)-subplans that satisfy Definition 5 are inserted into STS m+1 . In order to determine sp m+1 min , all the subplans generated by STS m should be scanned after they are generated. However, maintaining all the child subplans of STS m is a runtime overhead. sp m+1 min can be incrementally identified by tracing each of m-subplans in STS m . When sp m min is changed, the (m + 1)subplans already inserted into STS m+1 are re-examined. The pseudo-code of these steps is shown in Fig. 3. The following example shows how k-EGA finds the optimal plan of a multiway join query.
Example 2) Suppose that there are the SQL-like query Q in example 1 and the cost statistics in Table 2. Let k be 1.  erated by sp 1 1 . Likewise sp 2 3 and sp 2 4 are generated by sp 1 2 . The effective ranges of the subplans are presented in Fig. 2 (b). As is shown in Fig. 2 (b) is {J 1 , J 4 , J 8 } and sp 3 2 is {J 2 , J 6 , J 9 }. μ(sp 3 1 ) is lower than μ(sp 3 2 ), so that the final query plan becomes sp 3 The standard deviation of a subplan presents how much the subplan fluctuates. The small standard deviation means there is little change in the cost of the subplan. On the other hand, the large standard deviation means the cost of the subplan fluctuated widely in the past, and it is difficult to determine the current cost of the plan. Also, if any subplans fluctuate widely, the following subplan is also likely to be done because the result of the parent can be the input of the child. Therefore, k-EGA traces these subplans more deeply and intensively and this can cause a more optimal plan. (By Definitions 4 and 5, the subplans which have a large standard deviation are likely to be included in promising subplans.) As mentioned above, the number of simultaneously traced subplans is flexibly controlled. Actually, the capability of the k-EGA is from the general greedy optimization to exhausted optimization. If the range variable k is set to 0, only one subplan is traced since the minimum and maximum values of effective ranges of all subplans become the same by Definition 4. Consequently, it produces the same plan as the general greedy optimization. In contrast, if the range variable k is sufficiently large, so that all the m-subplans can be included in STS m , its resulting plan is the completely optimized plan.

Experiments
In this section, we provide the experimental results of the proposed algorithm. All the algorithms are implemented in C, and all the experiments are executed on a Pentium 4 CPU 2.66 GHz system with 1 GB RAM. The system runs Linux with a 2.4.5 kernel and gcc 3.3.2. To validate the cost model proposed in this paper, a multiway join query with 2 join operations, A B and B C, is experimented on for a synthetic dataset containing tuples of integer attributes. The domain size of each join attribute is 100. The value of each attribute ranges over [0..99] and is randomly selected. The input rates of the operand data streams A, B and C are 10 (tuples/s), 20 (tuples/s) and 50 (tuples/s), respectively. The cost of the query is compared with the actual execution time of the query by varying either the input rate of the operand A and the domain size of the join attribute of A B. Figure 4 shows that the cost is closely proportional to the actual execution time.
To analyze the performance of the k-EGA, a synthetic dataset is generated. For a multiway join query, the sizes of a sliding window and a hash table are fixed. The rates of the operand data streams are uniformly chosen from 10 tuples/s to 1000 tuples/s and the domain size of the attribute of each join operation is also randomly varied from 1 to 1000. The value of each join attribute is randomly selected from 1 to the domain size of the join attribute. Figure 5 shows the effects of a range variable k as the number of operand data streams for a multiway join query is increased. If the range variable k is set to 0, k-EGA is the same as the routing policy for a STAIR [7], which routes an incoming tuple to the join operation with the lowest selectivity. Intuitively, the execution time should increases as the number of join operations in a query increases. However, as k is increased, the execution time decreases since the final plan produced by the k-EGA becomes more efficient. In Figs. 6 and 7, a multiway join query with six join operations is experimented on for the nine synthetic datasets in Table 3 to compare the performance of k-EGA. The input rate and domain size of an operand data stream are randomly selected within the predefined range. The join selectivity of a join operation is determined by the reciprocal of the domain size of the attribute in each join predicate. In Table 3, the term a±b means that a value is randomly generated from the range [a − b, a + b]. For example, the input rates of the cost in DS 5 are randomly generated from the range [50,150] and the domain size of each join attribute is randomly generated from the range [500, 1500].
Given the k for a multiway join query, let Cost k denote the cost of the execution plan generated by k-EGA and Cost opt denote the cost of the optimal plan. To observe the performance of k-EGA with respect to the range variable k, we define the relative effectiveness E(k) of k-EGA. It is defined as.
The smaller the E(k) is, the closer to the optimal plan the execution plan generated by k-EGA is. In Fig. 6, E(k) decreases as k increases even if the optimality of k-EGA decreases as the range of the input rates becomes wider. As shown in Fig. 7, E(k) decreases as k logarithmically increases even if the optimality of our algorithm decreases as the range of the domain size increases.
To reoptimize the current execution plan continuously, it is important to monitor its cost periodically. If the plan is no longer efficient, k-EGA generates a newly optimized plan. Given the current plan p i , let Cost lastest (p i ) denote the cost of p i when the plan was originally generated. On the other hand, let Cost cur (p i ) denote the cost of p i at present. A newly optimized plan should be generated by k-EGA when the following equation is satisfied.
A cost allowable ratio δ is a user-defined parameter that denotes the allowable ratio of the current cost relatively to the latest optimized cost. As δ is set to be smaller, k-EGA changes its optimized plan more sensitively. Figure 8 shows how k-EGA adaptively reoptimizes its     execution plan for a seven-way join query. To simulate the dynamic change in the characteristics of input data streams, the tuples of the three datasets D 1 , D 2 and D 3 in Table 4   Table 4 Sub-datasets for Fig. 8.
are processed in sequence. The y-axis in Fig. 8 shows the cost per 1000 input tuples. The term "Optimal for D 1 " denotes the static optimal plan for the dataset D 1 only. When the tuples in D 2 start being processed, the cost of the optimal plan for D 1 sharply increases since the current evaluation sequence is no longer optimal to D 2 . However, k-EGA (k = 1) reoptimizes to generate a newly adjusted plan, and its execution plan cost remains to be close to the cost of executing the optimal plan. Figure 9 shows the trade-off between the optimality of the final plan generated by k-EGA and the time to generate the final plan. When k < 2, the cost of the execution plan optimized by k-EGA markedly decreases while the optimization time almost linearly increases. However, when k is large, the time to generate the final plan exponentially increases while the cost of the plan almost does not decrease. This means that when k becomes larger than a specific threshold (in this experiment, it is 2), the overhead of tracing multiple subplans simultaneously becomes greater than the profits gain obtained by finding a better plan.

Conclusions
In this paper, to optimize a multiway join query in a data stream environment, the K-Extended Greedy Algorithm (k-EGA) is proposed. The proposed algorithm traces a set of promising subplans simultaneously to reduce the possibility of producing a suboptimal plan. We provide a method of controlling the number of simultaneously traced subplans flexibly with respect to the range variable k. As a result, it is possible to adaptively control the number of subplans by changing the range variable k depending on the current processing load. Experiment results show how efficiently k-EGA can optimize a multiway join query in various stream environments. k-EGA can provide a good trade-off between the optimality of a generated execution plan and its optimization time. In addition, the plan generated by an MJoin operator cannot be utilized in a multiple query operation since the join orders of all operand streams are different. However, it is easy to extract common join operations in a number of plans generated by k-EGA. Most of the continuous queries are registered in advanced and should be processed continuously, so multiple query optimization is very inportant in a data stream environment.