{'fraction_non_alphanumeric': 0.03329788753630549, 'fraction_numerical': 0.010982906661164632, 'mean_word_length': 4.652462545976496, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Processing Pipeline for High Volume Pulsar Data Streams

Pulsar data analysis pipelines have historically been comprised of bespoke software systems, supporting the off-line analysis of data. However modern data acquisition systems are making off-line analyses impractical. They often output multiple simultaneous high volume data streams, significantly increasing data capture rates. This leads to the accumulation of large data volumes, which are prohibitively expensive to retain. To maintain processing capabilities when off-line analysis becomes infeasible due to cost, requires a shift to on-line data processing. This paper makes four contributions facilitating this shift with respect to the search for radio pulsars: i) it characterises for the modern era, the key components of a pulsar search science (not signal processing) pipeline, ii) it examines the feasibility of implementing on-line pulsar search via existing tools, iii) problems preventing an easy transition to on-line search are identified and explained, and finally iv) it provides the design for a new prototype pipeline capable of overcoming such problems. Realised using Commercial off-the-shelf (COTS) software components, the deployable system is open source, simple, scalable, and cheap to produce. It has the potential to achieve pulsar search design requirements for the Square Kilometre Array (SKA), illustrated via testing under simulated SKA loads.


Introduction
State-of-the-art pulsar search pipelines are comprised of two principal components. The first is a signal processing system. It converts the voltages induced in the receiving element of a radio telescope into digital signals, and identifies those 'significant' detections rising above the noise background. The signal processing system corrects for phenomena such as signal dispersion, excises radio frequency interference (RFI), and optimises search parameters yielding higher signal-to-noise ratio (S/N) detections. The signal processor ultimately produces some number of pulsar 'candidates' for analysis. These are time and frequency averaged data products describing each detection. For a complete description of the search process refer to [56].
The second search component is a filtering system. It identifies those candidates most likely arising from legitimate astrophysical phenomena, background noise, or terrestrial RFI. In principal this system allows signals of legitimate scientific interest to be isolated and set aside. For clarity we refer to the first component as the signal processor (SP), and the second as the data processor (DP). This paper focuses on the development of a new data processor. In particular we contribute the first realisation of a Commercial off-the-shelf (COTS) based DP for pulsar data designed to operate within an incremental data 79,12,96]. Emerging tools are often tailored to achieve specific science goals (e.g. improve sensitivity to longer period pulsars), and are often written with specific observing set-ups, or instruments in mind [issues reviewed and described in 61,62]. There are multiple DP systems in use at any one time, spanning global pulsar search efforts.
Almost all existing DP back-ends execute tasks sequentially, on batches of observational data [61,62]. The batches are usually processed off-line [77], either upon completion of a survey, or at the end of an observational session [e.g. 2,81,23,15]. A small number of real-time pipelines have appeared. These are designed to detect transient events [66,83], fast radio bursts (FRBs) [57,44,51,20,74,43], pulsars [69], or a variety of phenomena [10,26,21]. The adoption of real-time pipelines has occurred due to changing science requirements (desire/necessity for rapid follow-up), and in response to data storage pressures. These pressures are characterised by increasing data capture rates, which in turn yield increasing volumes of data [58,59,62,77]. Data pressures make it increasingly difficult to process data off-line, primarily as it becomes too costly to buy the media needed to store the data [77,62]. Such challenges are now impacting the search for pulsars. Without the capacity to continually store new data, search methods must be redesigned to maintain their effectiveness as the transition to real-time processing occurs [9,27,88,26].

Our Contributions
This paper makes four contributions in this area. Sections 2 and 3 introduce the core research problem, and describe the standard components of an off-line search pipeline (first contribution). After considering the feasibility of transitioning to on-line pulsar search (second contribution), and setting out the associated problems (third contribution) Sections 4-5 consider software tools useful for building a new on-line pipeline. Section 6 goes further, and presents a prototype pipeline design (final contribution). Sections 7-9 evaluate the prototype, and present results describing its performance. Finally Section 10 concludes the paper, and reviews the work.

Off-line Problem Definition
The goal for a DP is to reduce a set of candidate detections C, to a subset of promising candidates C ′ . These are the detections expected to possess the most scientific utility. The input set contains N elements, and N varies according to the exact processing configuration used 1 . The DP should return a set C ′ , such that |C ′ | ≤ |C|, though we desire |C ′ | ≪ |C|. 1 For one observation N ≥ 1000, for a survey N > 10 6 is normal.
Each element in C is describable as a candidate tuple. A tuple is a list of m elements. An individual tuple is defined as c i = {c 1 i , . . . , c m i }. Here each element is uniquely identifiable in C via the index i. For all c i ∈ C, it holds that |c i | > 0. For simplicity all c j i ∈ R, and there is no implicit ordering in C. In practice the numerical components of a tuple represent a detection, or its derived characteristics. The standard components usually stored within c i , are shown in Figure 1. The integrated pulse profile shown in a), is an array of continuous variables describing a longitude-resolved version of the signal averaged in time and frequency. The DM curve shows the relationship between candidate S/N and the dispersion measure (DM) [56]. Persistence of the detection throughout the time and frequency domains is shown in c) and d).
Persistence in frequency is represented by a two-dimensional matrix c), showing pulse profiles integrated in time for a set of averaged frequency channels (i.e. not full frequency resolution). Persistence through time is represented by a two-dimensional matrix d), showing the pulse profile integrated across similarly averaged frequency channels as a function of time. Other characteristics are also recorded. These can include the S/N, the DM, pulse period, pulse width, beam number, and acceleration; though other metrics are used [61]. The DP uses these characteristics to make accurate filtering decisions, producing as pure an output set as possible.
We note that C is not a multi-set, i.e., it does not contain exact duplicates. There are however non-exact duplicates present. These include pulsars detected by different telescope beams, or at slightly different characteristic values. Only the 'strongest' detection need be retained, with sub-optimal detections discarded to minimise |C ′ |. The strongest detection usually has the highest S/N. When C is available off-line, the strongest detection is found via exhaustive comparisons in the worst case.

Existing Data Processors
Existing pipelines execute filtering tasks sequentially upon C. The data is processed either after the completion of a pulsar survey, or at the end of an individual observational session if a faster pace of discovery is desired [see for example, 11,73,23,15]. At present it is feasible to permanently store all candidates in C. This allows detections to be processed more than once. It is common for new pulsars to be found in C, even after it has been searched multiple times [e.g. 45,67]. This happens when improved search algorithms/parameters are applied, revealing previously hidden detections.
There are some components common to all pulsar DP pipelines. These include 'sifting', known source matching, feature extraction, candidate classification and candidate selection. Together these successively filter candidates producing a purer set for analysis. The key components are  described in the sections that follow, so their function and computational requirements are understood.

Sifting
Sifting is an off-line matching problem. The goal is to to accurately identify duplicate detections in a candidate data set (i.e. find harmonically related duplicates), via comparing candidate pairs. For each arbitrary pair c i and c k , where i = k, a naïve sift will exhaustively compare all possible pairs using a similarity measure s. The measure is normally associated with a decision threshold t, applied over one or more variables in a candidate tuple. If s is above some threshold, the pairing is considered a match. Otherwise the pair is considered disjoint.
The accuracy of the similarity measure can be quantified, when the ground truth matching is known. In such cases the performance of s can be measured using a simple metric p, such as accuracy of the output matching, The current 'best' sift implementation [e.g. in the Seek tool, 55] performs an optimised comparison of each c i , to every c k in C. This has a memory complexity of O(n) as all n candidates must be stored in memory. The corresponding run time is approximately O(n 2 ). We note that minor adjustments to the approach can yield better runtime performance. Using combinatorics we find that a decreasing number of comparisons only need be done for each c i . The total number of permutations for a set of length n, where k items are compared at a time, is given by the binomial coefficient, This approach is dominated by O(n!) for all k. In practice for the required k = 2 the runtime is dominated by O( 1 2 (n − 1)n), an improvement over the worst case.

Feature Extraction
Numerical variables known as 'features' are usually extracted from candidates post-sifting. These are useful for deriving accurate filtering decisions, and are stored within each candidate tuple. Historically, features were comprised of standard signal characteristics (DM, pulse width, etc) used by human experts to filter candidates manually. Today features are utilised principally by ML algorithms [see 32,16,76]. These build mathematical models able to filter and separate candidate data automatically, based on their feature values. For pipelines employing ML tools, features are generally more complicated than simple signal characteristics (see Section 2.1.3 for more details). The computational cost of extracting features varies. The costs are rarely reported, though can be very high [61].

Candidate Selection
Candidate selection applies filtering decisions using the features that comprise a candidate tuple. Candidates passing through this stage are 'selected' for further processing or study. Candidate selection varies in complexity. From a single threshold applied over S/Ns, to far more complex automated machine learning-based approaches [e.g 33,13,97,68,62]. Where only thresholds are applied, runtime and memory complexity is constant per candidate. For machine learning and other more sophisticated methods, complexities are hard to determine as they are input data dependent. Candidate selection is notoriously difficult. In recent years a rise in candidate volumes has spawned what has become known as the 'candidate selection problem' [33].

Known Source Matching
Known source matching involves determining which detections correspond to known pulsar sources. The procedure is often carried out during sifting, though can be done independently after candidate selection. It requires a set K of known sources, usually obtained from a known pulsar catalogue [such as , 65]. For simplicity 2 each source in K is defined as k i = {k 1 i , . . . , k m i }, with each tuple uniquely identifiable in K via the index i. As before for candidates, all k j i ∈ R, with the meaning of each k j i the same as for candidates (or at least mappable). A brute force matching approach compares each candidate c i to every k i ∈ K. This corresponds to a runtime complexity of O(n·|K|). As new pulsars (and other possible radio sources) continue to be found over time, |K| is gradually increasing. The brute force method is thus computationally expensive 3 for an increasing |K| and an unbounded n. The memory complexity is O(|K|), as each known source is typically stored in memory to facilitate fast matching.

Manual Analysis
Following the previous steps, a set of candidates C ′ will be stored for manual analysis. Here experts manually examine the elements of C ′ , and judge their discovery potential. Given the large number of candidates usually in C ′ , this process is time consuming. This in turn introduces the possibility for human error [62]. Manual processing steps may have to be re-run when optimised search parameters/improved search methods are applied to C ′ .

Feasibility of Real-time Search
Empirical data describe a trend for increasing survey data capture rates over time [62]. This is expected to continue in to the Square Kilometre Array (SKA) era [80,30,17,61]. The projected data rates are large enough to make the storage of all raw observational data impossible (e.g. filterbank/voltage data), and the storage of all reduced data products (e.g. time vs. frequency vs. phase data cubes) impractical [62]. Without the capacity to store data for off-line analysis, SKA pulsar searches will have to be done in real-time at SKA-scales [30,77,61].
To overcome this problem, selection methods must be adapted to the real-time paradigm. A real-time DP must operate within the time and resource constraints imposed by instruments such as the SKA [28,29,72,82]. DP operations must also execute functionally within those constraints. This is difficult to achieve. Not all DP operations, as currently deployed, can run within a real-time environment. The following problems reduce the feasibility of transitioning off-line pipelines to the on-line paradigm: • Sifting is an inherently off-line operation that requires data to be aggregated. This is problematic in a real-time environment. The time taken to aggregate data can violate latency requirements. Whilst the memory required to retain all N candidates may not be available. For observations where rapid followup is crucial (fast transients), the delay induced via aggregation delays a rapid response.
• Pulsar feature extraction has never been undertaken on-line, and is currently an off-line batch process.
• Known source matching, similar to sifting, is currently executed over all candidates off-line. It therefore suffers from similar issues.
• The most accurate selection systems in use today, are built using off-line machine learning algorithms. Such off-line systems are inherently insensitive to distributional changes in the input data over time. Yet sensitivity to change is an important characteristic for an on-line system to possess, as it enables automated learning and operation.
Given such challenges, existing methods must be redesigned, or new techniques developed, to maintain the effectiveness of DPs for future pulsar searches. In some cases converting pipeline tasks to run on-line is deceptively trivial. For instance sifting could be done on-line, via only comparing each c i to c i−1 and c i+1 . Whilst functionally possible, such an approach would significantly reduce sifting accuracy. This would in turn produce many more false positive detections, increase the lead-time to discovery, and possibly preclude some discoveries being made.

On-line Problem Definition
In a data streaming environment, only a portion of C is available for processing at any given time. Either due to real-time constraints forcing batches of data to be processed separately, or C being too large to process in a single pass. In either case there are two processing models that can be employed [discussed elsewhere e.g., 98,19]. The incremental 'tuple-at-a-time' model applies when individual data items arrive at discrete time steps. An item c i arriving at time i, is always processed after an item arriving at time i − 1, and before i + 1. The batch (or microbatch/discretized streams) model applies when groups of items arrive together at discrete time steps. Here batch B i arrives after B i−1 and before B i+1 . This is summarised in Figure 2 for clarity. Both models temporally order data, which has implications for how the data can be processed.
The batch model is advantageous when groups of items exhibit distributional similarity. For instance when processing data from seismic sensors, it makes sense to work with batches, as patterns of seismic activity will likely be close temporally. Whilst if looking for fraudulent activity in a stream of random financial transactions, each transaction should be processed in isolation. Otherwise an innocent transaction may be incorrectly linked to a fraudulent one. There is a trade-off between computational efficiency and practical utility [22,54], that must be struck when choosing a processing model.
For the processing of single pulse events, it does not make sense to batch data. Waiting for a batch incurs a time penalty which impedes our ability to initiate rapid follow-up. For pulsar search batch processing is possible, however batch sizes are potentially very large. In both cases it would appear that an incremental model is the simplest to adopt.

An Incremental Model for Data Processing
We extend the notation used previously. The candidate set C is now a candidate stream C = {c 1 , c 2 , . . . , c n , . . . }.
Here each c i is delivered at time step i, and c i always arrives before c i+1 (discrete time model with temporal ordering). The size of the set C is unbounded in this scenario (unknown N ). This assumes a limitless supply of individual candidates, redefining the SP as a data stream producer, and the DP as a data stream consumer.

Practical Assumptions
We assume a worst case scenario, where each data item is processed just once, as per the incremental model. Filtering decisions must be made on an individual candidate basis, with limited knowledge. For example, suppose c 2 arrives for processing where only c 1 has been observed. Here a filtering decision can only be made using knowledge of the candidate/s already seen and c 2 . If c 2 is most likely noise or interference, then the decision is conceptually simple. Yet the decision making process becomes complicated in scenarios where multiple detections of the same source are likely to be made during an observation. If c 2 is indeed a pulsar detection, do we retain it? Is it the strongest detection we are likely to see? If it is weak, should we throw it away, assuming a stronger detection will come along? Do we keep every weak detection, and risk increasing the size of C ′ ? There are no simple answers to such questions.
To proceed in our domain, we assume candidates will be ordered according to c 0 i (an ordering variable). This ordering must be strict, so that ∀c i ∈ C, c 0 i ≤ c 0 i+1 . This definition specifies an ascending order, though it would make no difference if it were descending. By using an ordering we can have confidence that similar items will be close together in the stream (close temporally). This simplifies our data processing, and will be used to develop an on-line sifting algorithm in Section 6.2. Note that we are implicitly assuming there exists an upstream processing component, capable of ordering data items correctly. We accept that it may not always be possible to assign an ordering conducive to improved data processing. However, where possible we should attempt to order our data in a meaningful way, as it greatly assists with the downstream processing. We now look for frameworks that can accommodate this processing model.

Candidate Frameworks for Prototyping
The design of a DP is driven by multiple considerations. We have focused upon utilising COTS software components to reduce cost. This is our primary design driver. However our choice of software components was also driven by the need to reduce computational overheads, whilst maximising scalability and design modularity.

Review & Chosen Framework
Signal processing systems often utilise accelerator hardware to enable real-time operation. This includes Graphics Processing Units (GPUs) [e.g. 64,25,30] and Field Programmable Gate Arrays (FPGAs) [e.g. 48,31,95,90,91]. In contrast most DP software is executed on general purpose computing resources, with some exceptions 4 . We therefore only consider DP execution upon standard Central Processing Units (CPUs). There are many frameworks/approaches that meet our criteria [for a review see e.g., 35]. Yet we chose Apache Storm [39,4] to support our new on-line prototype. It was chosen due to its application in the real-world to similar large-scale incremental processing problems [see 41,86,92] and strong user support base. Storm is similar to frameworks such as Apache Hadoop [94,40,7], Samza [49,6], Spark [99,8] and S4 [71,5], and is fairly representative of similar modern COTS tools. Since completing this work, we discovered other frameworks in the software ecosystem similar to Storm. These include Apache Samoa [14,3], Apache Flink [19], and Kafka Streams [78]. The use of Apache Storm in lieu of these systems, does not undermine the novelty or usefulness of the contributions of this work. Principally as the processing pipeline proposed in Section 6., is framework independent and compatible with any tool capable of supporting the 'tuple-at-a-time' processing model.

Apache Storm-based Processing Framework
Storm is a Java-based framework. It operates under an incremental data stream model. It is underpinned by the notion of a directed-acyclic graph [DAG, see 85]. The graph models both the processing steps to be completed, and the data flow. The graph is known within Storm as a topology. Storm topologies only allow data to flow in one direction. This makes recursive/reciprocal processing steps impractical to implement (though newer frameworks discussed in Section 4 can overcome such limitations).

Topologies
A topology is comprised of nodes. Nodes represent either data output sources known as spouts, or processing tasks applied to the data called bolts. Data flows from spouts to bolts via edges in the topological graph. Individual data items are transmitted in the form of n-tuples. These are finite ordered lists much like the candidate tuples defined in Section 2.

Edges
Edges represent generic connections between processing units, made via a local area network (LAN), or a wide area network (WAN). Data flows via the edges as tuples. Upon arriving at a bolt a tuple is usually modified, and the 4 Machine learning based filters [e.g. 97]. updated tuple emitted for a downstream bolt to process. The practical data rates between edges are non-uniform, and often subject to change. Fluctuating data rates can lead to resource contention, if too much data is funnelled through too few processing bolts.

Bolts
Each bolt encapsulates a modular computing task, completable without interaction with other bolts. Such modularity allows data to be processed in parallel, across heterogeneous computing resources without a loss in functionality. The duplication of bolts allows the topology to be scaled to increasing amounts of data, without additional development effort. It also means that when a bolt fails, it can quickly be replaced with a new instance. The modular design is not ideal for all scenarios. Particularly those where the persistence of state is required. Generally state persists within a spout/bolt only whilst it is executing. If it fails or is stopped, that state is lost. Storm more generally is unsuitable for processing components that are i) not lightweight, ii) require a persisted state to function correctly, or iii) require data to be aggregated/batched.

Tuple Flow
It is possible to control the flow of tuples through a topology. This is achieved via a cardinality relation. This is defined between spouts and bolts via edges. The relation can be one-to-one, one-to-many, or one-to-all. If using a one-to-one relation, tuples either move forward to a single random bolt, or a value-specific bolt based on a tuple attribute test 5 . Upon transmission of a tuple, a 'tuple received' acknowledgement can be requested by the sender (originating spout or bolt). If the sender receives no acknowledgement, it will resend the tuple. This enables Storm to maintain the property that all data is processed at least once, helping ensure fault tolerance.

Deployment
Topologies execute upon clusters comprised of two distinct nodes. There are i) master nodes that run a daemon called 'Nimbus', and ii) one or more worker nodes running a daemon called 'Supervisor'. A master node is responsible for executing tasks on the worker nodes, and restarting spouts/bolts after failures. The worker nodes execute spouts and bolts separately, within threads spawned by a Java Virtual Machine 6 (JVM).
Communication between the master and worker nodes is coordinated by Apache Zookeeper [40]. Zookeeper maintains any state required by the master and the supervisors, and restores that state in the event of failure. A Storm cluster therefore consists of three components: Nimbus, one or more Supervisors, and a Zookeeper. The prototype DP topology. Data is supplied to the topology via SP spouts, which generate candidate tuples. The tuples are propagated forwards through the topology. No tuple is duplicated, and each reaches only one bolt in the topology at any given time. There are nine distinct types of bolt in the topology, described in the sections that follow.

DP Prototype Design
The prototype design is shown in Figure 3. It has a single 'layer' of input spouts, which emit candidate tuples at a controllable rate. Candidates emitted from the spouts are propagated to a single random bolt, at each subsequent layer in the topology. The first layer contains data ingest bolts. These modify the data so that it is transmitted in a format amenable to further processing.
The second layer attempts to filter out duplicate candidates via sifting. To achieve this, the sifting bolts employ a new distributed global sift algorithm. This is described in more detail in Section 6.2. After sifting, tuples are sent to bolts that perform data pre-processing (i.e. normalisation). This is required prior to ML feature extraction. Once complete, the next layer of bolts extracts the features. The features are appended to the tuple, and sent onwards for further processing.
Following feature extraction, tuples are passed to ML classification bolts. These execute an on-line ML algorithm developed by [62]. This predicts the true class origin of each candidate. The predictions are appended to each tuple, and passed on to secondary sifting bolts. These remove duplicates in lieu of the additional information obtained during ML classification. Tuples making it through the second sift, are passed to source matching bolts. These attempt to match promising candidates to known pulsar sources in a pulsar catalogue. Candidates not matched to known sources (likely new pulsars) generate alerts for follow-up action.
The topology has been designed to provide auditing capabilities. Auditing nodes are connected to the processing bolts, via the dashed lines shown in Figure 3. They are used to audit the performance of individual bolts, or the entire topology during testing. The auditing nodes are not activated during normal execution.
The individual spouts/bolts are now described in more detail. Source code for the procedures described in the following sections can be found on-line [63]. The code is provided in the form of an interactive iPython notebook. It contains additional details that support the information presented in this paper.

SP Spouts
At present no SP system can generate data fast enough to stress our prototype. We therefore created SP spouts which emit real candidate tuples driving the data processing. These play the role of the SP, delivering data at a controllable rate. The spouts are highly customizable. They can generate a fixed number of candidates, or as many as possible within a fixed period of time. Importantly the ground truth label for each and every generated tuple is retained. This allows candidate filtering accuracy to be evaluated anywhere within the topology.
To be realistic, the spouts must output real-world candidate class distributions. For the prototype we define a ratio used to achieve this. Suppose each candidate c i is associated with a label defining its true origin. This can be modelled numerically via a binary variable y, where describes the imbalance between pulsar and non-pulsar candidates in generated data. The ratio observed during real pulsar searches varies from 1:7,500 [46,84] to 1:33,000 [24,52,1]. For the prototype, we maintain ratios of up to 1:10,000 (0.0001). Evidence suggests this figure to be representative [62] and challenging to deal with.
To maintain the ratio, spouts have been designed to emit two types of candidate. Type 1 candidates are contrived randomly generated non-pulsar tuples. Whilst valid, they simply resemble white noise. Type 2 candidates are representations of real candidates. These were sampled from data obtained during the High Time Resolution Universe Survey South [HTRU, 46,60]. Type 2 candidates provide a genuine test of our prototype's discriminative capabilities. Type 2 candidates can describe either the pulsar or non-pulsar class. As the sample of type 2 candidates is small, some pass through the pipeline more than once during testing. This does not invalidate our results, as duplicates must be passed through the system to test our sifting approach.

Spout Data Rate and Volume
Input data rates are determined by the length of an observation t obs (seconds), the number of telescope beams n beam , the number of candidates returned by the SP per beam c beam , and the size of each candidate c size (MB). The total data volume per observation is given by, Whilst the SP to DP data rate, assuming a steady uniform transmission, is given by, The assumed parameter values are given in Table 1. These were chosen after studying SKA design documentation [see 28,29,17]. Here candidate tuples are comprised of, • a 262,144 sample data cube (128 bins, 64 channels, 32 sub-ints), where each sample is 8 bytes in size. This describes the detection in time, phase, and frequency.
Together these components produce a tuple ≈ 2.2 MB in size. For 1,500 beams, and 1,000 candidates per beam, this equates to 1.5 million candidates. The corresponding total data volume is ≈ 3.3 TB per observation. The t obs = 600 second time constraint, implies a data rate of 5.5 GB/s. The SP spouts will generate data at and above this rate, when simulating pulsar search operations.

On-line Global Sifting
Standard sift only compares candidates detected within the same telescope beam [e.g. 42]. Our approach works globally across beams, in principle yielding better results. It sifts one candidate at a time, as opposed to candidate batches. It's success is predicated on two assumptions. First, it assumes a stream partitioning that ensures candidates with similar periods from all beams (measured in µs) arrive at the same bolts. Thus an individual bolt represents a period range, and together they cover the plausible range of pulse periods. Second, it assumes candidates are ordered (see Section 3.2) according to their pulse periods. This facilitates the partitioning. Since similar period candidates always arrive at the same bolts, it becomes possible to check for duplicates via counting observed periods. The logic underpinning this approach is simple. If a period is observed many times, it is a possible duplicate. When possible duplicates are compared via other variables (e.g. DM) for additional rigour, counting can be used to efficiently find duplicates.
The general approach is summarised in Algorithm 1. This describes the code at a single bolt. Note the algorithm requires sufficient bins to count accurately. As pulsar periods are known to be as low as 1.396 milliseconds [36], and given that pulse periods can be similar at the millisecond level; there must be enough bins to count at microsecond resolution. The algorithm uses the array F to maintain the count, by mapping the array indexes 0, . . . , n to specific periods/period ranges. We use a scaling and a rounding operation to achieve this. Together these find the correct bin index to increment for an arbitrary period. The variables f loor and ceiling help to accomplish this. The scaling of the data is done on line 8, with the integer rounding operation done on line 9 which obtains the index.

Algorithm 1 Global Sift
Require: An input stream C = {..., (ci), ...}, such that each ci is a candidate, and c j i its j-th feature. Here c 0 i is the pulse period (ordering variable). Requires a similarity function s which can be user defined, the number of period bins p b to use, the smallest period value expected at the bolt min, and the largest period value expected at the bolt max. n ← n + 1 ⊲ increment observed count The if-else statement on lines 10-15 contains the only logic in the algorithm. It checks if the period counter is greater than zero (indicating it has already been observed). If F [index] > 0, then we pass this information to the similarity function for checking. The algorithm can be modified for improved accuracy. For example, if frequency counts are maintained for multiple variables (DM, pulse width, beam etc.), the probability of a match can be estimated across all of them. Note the similarity function s used on lines 12 and 15, can also be made to check matches counted in neighbouring bins. Candidates counted in bins neighbouring F [index] have very similar periods, and thus should be considered possible duplicates. It is trivial to implement, though requires additional memory and a sliding window of examples. The use of a sliding window is discussed briefly below.
The general approach described has weaknesses. The first candidates arriving at a sift bolt will never be flagged as duplicates, due to no prior periods being been observed. Furthermore, the approach will treat duplicates with different S/Ns as equivalent. However we wish to retain only the highest S/N detection, and discard the rest. The model described thus far cannot achieve this. It cannot anticipate if, or when, a higher S/N version of a candidate will enter the topology. Such weaknesses are a result of the tradeoff between computational efficiency, and sifting accuracy. This problem can however be overcome, using a sliding 'similarity' window over the data [see 63].

On-line Feature Extraction
Eight features are extracted from each tuple moving through the topology. These are described in Table 4. of [62]. The first four are statistics obtained from the integrated pulse profile (folded profile) shown in plot a) of Figure 1. The remaining four similarly obtained from the DM-SNR curve shown in plot b) of Figure 1.
The computational cost of generating the features in Floating-point operations (FLOP) is extremely low. If n represents the number of bins in the integrated profile and DM curve, the cost is as follows: • Mean cost: 2(n-1) + 4 FLOP.
This assumes the following: addition, subtraction and multiplication all require 1 FLOP, whilst division and square root calculations require 4. The Standard deviation calculation assumes the mean has already been calculated first. Likewise, the skew and kurtosis calculations assume the mean and standard deviations have already been computed, and are simply reused. The features cost 2,848 FLOP per candidate [61]. The runtime complexity of the feature generation code is O(n) over the input space, and memory complexity is O(n). Note that in both cases n is small.

On-line Classification
This prototype is focused on processing incremental data streams in real-time. At present only one classifier has been developed specifically for incremental pulsar data streams -the GH-VFDT [59,62,61]. All other pulsar classifiers developed in recent years can be characterised as off-line supervised systems. Thus in terms of classifiers, the GH-VFDT is the logical choice for this work. Despite this we acknowledge that we may have overlooked classifiers capable of achieving better recall on our data. As the goal of this paper is to illustrate the feasibility of processing pulsar data incrementally in real-time, and not to advocate the use of a specific ML algorithm, this is not something we consider to be problematic. Our design does not preclude the use of any other incremental data stream classifier developed now or in the future.
The GH-VFDT classifier was designed for pulsar candidate selection over SKA-scale data streams. It is an on-line algorithm, capable of learning incrementally over time. It is therefore able to adapt to changing data distributions, and incorporate new information as it becomes available. The GH-VFDT is extremely runtime efficient, as it was designed to utilise minimal computational resources.
The algorithm employs tree learning to classify candidates, described via their features [for details of ML classification see, 32,16,76]. Given some ground truth 'training data' describing the pulsar and non-pulsar classes, tree learning partitions the data using feature split-point tests [see Figure 8 in 62]. Split points are chosen that maximise the separation between the classes. In practice this involves firstly choosing a variable that acts as the best class separator. Once such a separator is found, a numerical threshold 'test-point' is found, that yields the greatest class separability. This process repeats recursively, producing a tree-like structure. The branches of the tree form decision paths. The paths can yield highly accurate classification decisions when given high quality training data.
The memory complexity of the algorithm is O(lf · 2c) (sub-linear in n), where l describes the number of nodes in the tree, f the number of candidate features used (f = 8 for the prototype), and finally c the number of classes (here c = 2, pulsar & non-pulsar). The runtime complexity of the algorithm is difficult to quantify, as it is input data dependent. However it is of the order O(n).

On-line Post Classification Sift
Post-classification sifting aims to remove duplicate detections in lieu of ML predicted class labels. Here two candidates with very similar period and DM values (other variables can be considered), have their predicted class labels compared. If the same, these are considered likely duplicates. In this case, only the strongest detection need be forwarded on. This approach can be improved if a sliding window over tuples is used [see 63]. This allows similar tuples to be compared together in small micro-batches defined by the window. Only after the window moves away from an unmatched tuple does it get forwarded.

On-line Known Source Matching
We have developed a new, very fast, source matching algorithm. It is built upon a divide and conquer programming approach (built upon tree-search). This involves attempting to recursively divide the matching space, reducing the number of comparisons to be undertaken on each iteration. It relies on an ordering applied over the set of known sources K, represented as an array. A total ordering of elements in K is required, according to some variable k j i . To achieve a total ordering, then for all k j i , k j i+1 , and Here equations 6-9 define the antisymmetry (6), transitivity (7), reflexive (8) and totality properties (9) respectively. To apply the desired ordering, we require a numerical value per source that satisfies these properties. This can be obtained via measuring the angular separation θ, between each known source, and a single reference coordinate (00 h 00 m 00 s and 00 • 00 ′ 00 ′′ ). This allows sources to be strictly ordered according to their separation from the reference point. This reference value should be computed off-line so that the array K is ordered correctly in advance. Note that here we assume K contains all known sources. This represents a worst case matching scenario. In practice, K will likely be filtered, so that it includes only those known sources in the patch of sky in which the candidates were observed (pre-selected known sources based on pointing location).
For each candidate source c i to be matched, the reference separation must also be computed. Intuitively known sources near to c i would appear to be similarly separated from the reference point. The reality is more nuanced. It is possible for sources to be similarly separated from the reference point, yet be very far apart. This is shown in Figure 7. Here candidate c i and known source k 1 , are both separated from the reference point by an angular separation of θ. However these are too distant from one another to be considered a match. This does not affect the matching accuracy of the new approach. The sorting of known sources is only used to reduce the search space. In particular it allows us to find a position from where to start searching in the known source array K. We call this position the 'search index'. From there, we need only compare c i to sources in K around the search index, i.e. with separations ≤ 2θ with respect to the reference point. These known sources must be adjacent to the search index given the properties of the ordering.
A divide and conquer algorithm firstly divides up the search space recursively, until a search index is found. This is achieved as follows: if the separation of c i from the reference point, is less than the separation between the reference point and the centre source in K, search the first half of K. Else, search only the second half of K. This process is repeated until the closest single match is found, and a search index returned. Once the closest match is found, a secondary algorithm compares c i to those sources nearest to it. Precisely how many comparisons are done to the left, and right of the search index, is up to the user. For the prototype, comparisons are not undertaken if the distance between c i and k i exceeds 1.5 • .
The divide and conquer approach is presented in Algorithm 2, and the matching procedure defined in Algorithm 3. The matching procedure compares the period and DM of a promising c i , to some potential match k i . The known source k i is considered a possible match for c i , only if their period and DM values are similar to within 3.
1. Figure 4: A visual representation of how the matching search space is determined. In a) we see candidate source c i , which has an angular separation from the reference point of θ. There are also 2 known sources k 1 and k 2 , with the same separation from the reference point. The angular separation between c i and other known sources is computed by Φ. If the separation is greater than θ, then c i is not compared to k i . Otherwise a comparison is performed. From this representation it is clear that c i could be related to known sources up to 2θ from the reference point. However only the region around c i should be searched. In b) we see the same information, but stored within an array. This is the array used for searching. The sources in the array are sorted according to their angular separation to the reference point. This ordering allows the search region to be found quickly in step 1. This is then refined, producing the narrowed search region in step 3.

Algorithm 2 Divide & Conquer Matching
Require: A set of known sources K = {ki, . . . , km}, the numerical sorting variable sep, of the candidate ci to compare. Finally istart is the index to search from in K, and i end is the index to search to. a user specified error margin e margin ∈ [0, 1]. For example, an e margin = 0.1 corresponds to a 10% error margin. When using this margin we consider a known source to be possible match, only if its period and DM are within 10% of the candidate's (±5%). A 10% error margin over pulse period is large. A smaller margin should be used in practice, especially for data recorded with high time resolution. However for our test data, pulsar periods differ by up to ±5% compared to values recorded in the pulsar catalogue. The prototype does allow for a different margin to be chosen.
The computational complexity of the overall approach is O(n · τ ), where τ is a proxy for the number of comparisons made between known sources and candidates based on θ. The modified runtime is practically speaking linear in n (as τ is usually small). The complete algorithm has been described in an iPython notebook [63]. The approach presented here is inherently simple and fast. It differs from existing cross matching methods used in astronomy, which employ Bayesian methods [18], or multiple catalogues describing sources at different wavelengths [34].

On-line Alert Generation
Prototype alert generation bolts have limited functionality. These do not generate genuine alerts, since alerts are not required for our prototyping efforts. Thus alert nodes  if c dm < k dm + dm dif f then 13: if c dm > k dm − dm dif f then 14: sep ← calcSep(ki, ci) 15: if sep < θ then 16: possibleM atch(ki, ci) simply act as processing units, that slow down computation as though alerts were being generated. Similarly as no archival system exists, archival is not simulated.

Auditing
The prototype incorporates processing nodes able to audit runtime/filtering performance. Auditing is accomplished in two ways. Per node auditing records the filtering accuracy and runtime performance, for a specific node only. This incurs a small computational overhead impacting the node being audited. As this type of auditing is only intended for use during development, it doesn't affect runtime performance during scale testing.
End-to-end auditing measures filtering accuracy and runtime performance across the pipeline. This is achieved without incurring additional runtime overheads via the use of time stamps. Unique time stamps are attached to each tuple upon entering and exiting the topology. The time stamps accompany tuples through the topology, and record time to microsecond resolution. By determining the difference between time stamps, an estimation of the time taken for an individual tuple to move through the topology can be determined. Individual tuple transfer times can also be aggregated. By averaging over all tuples reaching the end of the topology, average tuple processing times can be computed. Additional metrics are also monitored at the auditing bolts.

Simulations
Two forms of simulation were undertaken to test the prototype. The first consisted of small-scale simulations executed on a single machine. These were useful for testing and debugging the topology design and processing code. The second involved a larger scale deployment of the topology to the Amazon Elastic Compute Cloud (EC2). The cloud simulations were intended to assess the scalability of the system, and determine ease of deployment. In both scenarios the goal was to recreate the delivery of data from the SP to DP, during a plausible pulsar search scenario.

Local 'Cluster' Mode
Simulations were undertaken on single a computer running OSX 10.9. It possessed a single 2.2 GHz Quad Core mobile Intel Core i7-2720QM Processor, with a peak theoretical performance of 70.4 GFLOPs [38]. It was equipped with 16 GB of DDR3 RAM, and two 512GB Crucial MX100 solid state drives. A Storm cluster (version 0.95) was deployed on this machine, and the DP topology run in local cluster mode. This enabled testing of the framework prior to a larger scale deployment.

Cloud Infrastructure
We were awarded compute time upon Amazon's cloud infrastructure via the SKAO-AWS AstroCompute grant programme 7 . This time was used to test the performance and behaviour of the prototype, when scaled beyond a single machine. Using the Amazon Web Services (AWS) console, we provisioned a number of EC2 instances 8 . The provisioned instances are described in Table 2. Note it is difficult to estimate the overall compute capacity possessed by these cloud resources. This is because EC2 instances are deployed on shared hardware, subject to load balancing policies and stress from other EC2 users. Amazon describes the compute capacity of its virtual instances in terms of EC2 Compute Units (ECUs). According to Amazon's documentation, a single ECU corresponds to a 1.0 -1.
The Xeon E7310 (1 ECU) is capable of a theoretical throughput of approximately 25.6 GFLOPs, according to both Equation 10 and Intel's own export specifications 9 . The 241 ECUs used during cloud experimentation, therefore correspond to an approximate computational capacity of 6.2 TFLOPs.

Runtime Performance
Topology performance is measured using the auditing nodes. These maintain statistics which are updated after each tuple is processed. Crucially, performance results differ according to the topology configuration used. The configuration describes the number of spouts and bolts in each layer, and the cardinality relationship between them. We adopt a colon delimited notation to describe the configuration, e.g. 1 1... * : 2 1...1 : 1. This corresponds to one input spout in layer 1, two processing bolts in layer 2, and a lone bolt in layer 3. The superscript defines the cardinality relation between the current layer and the next. For this example, there is a one-to-many relation between the spout and the bolts in layer 2, and a many-to-one relation between the bolts in layer 2, and the lone bolt in layer 3.

Filtering Accuracy
The spouts ensure ground truth class labels are known for all candidates entering the topology a priori. It is therefore possible to evaluate filtering decisions at any bolt or spout. There are four outcomes for a binary filtering decision, where pulsars are considered positive (+), and non-pulsars negative (-). A true negative/positive, is a negative/positive candidate correctly filtered. Whilst a false negative/positive, is a negative/positive candidate incorrectly filtered. It is desirable to have as few false negative/positive outcomes as possible. The outcomes are eventually evaluated using standard metrics such as recall, precision, accuracy and F1 score. These are borrowed from machine learning research 10 [metrics used listed in 62]. Note that we do not compare classification results obtained with the GH-VFDT against existing off-line pulsar classifiers. This is because it is inappropriate to compare off-line and on-line methods in this way -they are intended 9 See http://www.intel.com/content/dam/support/us/en/documents/processors/xeon/sb/xeon_7300.pdf. 10 We do not use imbalanced metrics here (e.g. the G-Mean or Mathews Correlation Coefficient) as we are not trying to show superior imbalanced class performance, but rather pipeline feasibility).
for fundamentally different processing paradigms.
For the evaluation we use imbalance ratios much higher than those typically observed in the real-world. We do this for a simple reason -if pulsar examples are rare in our test data streams, the prototype will only experience minimal computational load following the classification step (i.e. most candidates will be filtered out here). Since known source matching is a relatively expensive computational operation, we wish to push the prototype to match many more candidates. Hence why we use c ratio = 0.05 or c ratio = 0.1, versus the real-world c ratio = 0.0001 (or worse).

Local 'Cluster'
The results of local cluster experimentation are given in Tables 3 and 5. The prototype is able to process 1.5 million candidates in well under 600 seconds, as shown in Table 3. It is able to achieve this at a very fast rate, exceeding 1,000 candidates processed per second (up to 6,000 per second during high throughput tests). Thus the prototype is capable of running functionally at SKA scales, on present day commodity hardware. It is worth noting this result was achieved on a single laptop, possessing only 70.4 GFLOPs [38] of computational capacity (theoretical max) 11 . This result appears valid, as pre-processing and feature generation (the most computationally intensive tasks) require only 19 MFLOPs of compute.
The prototype worked well functionally in local mode. The filtering accuracy was very high. It successfully removed most uninteresting candidates and duplicate detections. A 99% end-to-end accuracy level is achieved using our new global sift algorithm, with a corresponding high level of pulsar recall (81%). At the bolt level, the performance of global sift is also good, but not exceptional. Mistakes made by global sift account for why the pulsar recall rate did not move beyond 81%. During testing global sift returned 133 out of 222 unique pulsars injected. For a trivial algorithm which can be easily modified, this result is still promising. To illustrate this, consider the accuracy achieved when using a random sifting approach. Here a 50:50 sift is employed, that propagates 50% of tuples chosen at random. The accuracy achieved using a random sift drops to 84% whilst recall tends to zero.
Overall filtering accuracy is good, though below that achieved by recent ML methods used for candidate filtering [see 62,61]. This result is unsurprising and excepted,    since these methods are working with off-line data (models built off-line using large training sets). Yet we have still advanced the state-of-the-art, as these existing methods cannot be applied within the on-line scenarios we are considering. In these scenarios there is no training data available upfront, and severe computational restrictions.
Whilst the prototype's recall rate is high, it is not high enough for use with the SKA. Applied as is, it would miss almost 20% of all pulsars entering the system. However with improvements to sifting and known source matching possible, and the application of more sophisticated ML methods, this result can be greatly improved upon.
The individual bolts in the topology appear to be runtime efficient. Each tuple required on average only a few milliseconds of processing time at each bolt. The runtime per tuple does increase as more data is processed as shown in Table 3. The same is true of the total runtime. There is a disparity between the average processing times shown in Table 3, and the bolt processing times in Table 5. The disparity is accounted for by the time taken to transmit each tuple through the network. The values in Table 3 include this transmission time, whilst the values in Table  5 account for only the processing. The difference serves as a reminder that transmission time plays a significant role in the total runtime of a Storm-like system. It also emphasises the importance of running tests outside of local cluster mode, which does not incur the communication overheads experienced in the real-world.

Cloud Infrastructure
Cloud infrastructure simulation results are presented in Tables 4 and 6.  Here c obs is the total number of candidates entering the bolt, and tavg(ms) the time taken on average to process a tuple. Experiments run using c obs = 50, 000, and c ratio = 0.05 (produces many duplicates for sifting).
according to the total time it takes to process a tuple on average. The results in Table 4 show an increase in processing time, when initially beginning to scale the topology (when using 1-4 workers). This is followed by a decrease in processing time when more workers are added (increasing the computational power available during execution). Note we do not show filtering accuracy results in these tables. This is because filtering performance does not change between local and remote mode -the bolts and their filters are unchanged. The only difference is that remote mode allows the topology to be scaled to utilise greater computational resources, speeding up processing times and increasing capacity.
The initial increase in processing time observed when beginning to scale the topology is explainable by the results shown in Table 6. These describe the performance of the processing layers, as they are scaled-up (more bolts added) to process the data. When a topology contains only a few instances of each type of bolt, layers reach their pro-cessing capacity quickly. This causes an overall increase in processing times due to resource contention. As more bolts are added to the topology, counter-intuitively, the contention is not necessarily reduced. Whilst some bolts will begin to cope with the load when more instances are added, others will not. This happens when resource contention is shifted to another location in the topology, increasing processing times there instead. Note that these results were achieved when running the random 50:50 sifting approach described in Section 9.1. The random sifting approach had to be used instead of our global sift, as global sift is efficient enough to prevent us from stressing the downstream bolts significantly. As random sift propagates far more tuples, it allows us to study the scalability of the topology under more extreme loads.
Contention arises when bolts undertake different proportions of the total computational workload. In the case of the pulsar search topology, known source matching is the most computationally expensive procedure (see Table  6). Adding more bolts to the lower layers of the pulsar search topology, increases the tuple throughput reaching known source matching bolts. Initially the matching bolts cannot meet this demand. Only when at least 8 worker nodes are available, with enough known source matching bolts, does the contention disappear and processing time decrease. This is an important observation. The computational demands of known source matching are surprising. Table 6 shows which bolts experienced most resource contention. This is indicated by the capacity value (CV), CV = (tuples executed × avg. execute latency) measurement time . (11) A value of 1.0 corresponds to a bolt at full capacity. Values greater than 1.0 indicate a bolt over capacity, and less than 1.0 under capacity. When only two bolts are present in each layer of the topology, all bolts are at or near capacity. As the topology is scaled, most bolts begin to become under utilised. The exceptions are the ML classification bolts, and the known source matching bolts.
As more bolt instances are added to each layer of the topology, execution time latency reduces, and total runtime decreases. In some cases the runtime decrease scaled linearly with the number of worker nodes used. However, this is not always the case. There are fluctuations in the results which make it difficult to discern a genuine trend. The overall results shown in Table 4 indicate an improvement in performance scaling sub-linearly with the total number of worker nodes. This impression is based on averaged results, and agrees well with empirical experience. It is therefore likely the most accurate indication of true system scalability. We do not report total runtime, as all tests completed well within the t obs = 600 seconds required in the pulsar domain. This suggest our topology can process data in in real-time.

Conclusions
A prototype data processing pipeline for pulsar search has been developed. It is capable of on-line and realtime operation. It employs a combination of resource efficient algorithms and optimised selection methods. Together these enable large numbers of pulsar candidates to be filtered very accurately, using limited computational resources.
The performance of the prototype was first assessed on a single commodity laptop. During testing it was able to process 1.5 million pulsar candidates 12 , in under 600 seconds. It is therefore functionally capable of processing data fast enough to meet SKA design requirements (exceeds a processing rate of 6,000 candidates per second). The prototype was also deployed to a cloud-based software infrastructure. Here the system was similarly able to process data at SKA scales, using modest computational resources (6.2 TFLOPs of processing power). The runtime performance of the prototype scaled sub-linearly with the number of worker nodes used to execute the processing. Better scaling is likely impeded by the overhead of inter-node communication (i.e. latency and bandwidth restrictions incurred due to network communication).
However, the prototype was designed to favour computational efficiency over filtering accuracy and pulsar recall. This was done to ensure feasible operation at SKA scales, assuming a worst case processing scenario. This trade-off reduced pulsar recall. Thus although filtering accuracy reached 99%, the corresponding pulsar recall rate ranged between 81-88%. This is below the level required for SKA use. There is room to reverse the efficiency-accuracy tradeoff, without significantly compromising runtime performance. It would be reasonable to double the resource use of the system, to achieve higher filtering accuracy and pulsar recall.
Finally we emphasise that it was not our intention to suggest we should use so few resources (i.e. a modest 6.2 TFLOPs easily surpassed by a single modern GPU 13 ) to run a pulsar search pipeline. Rather we aimed to show that we can do so, using COTS tools at a low price point. A recall rate of 81-88% is promising, but the remaining 12-20% is incredibly difficult to isolate. As the recall rate increases, it becomes increasingly difficult to improve it further. Methods capable of isolating the last few percent will likely be sophisticated, and possibly be accompanied by higher computational runtime costs.
We recommend that future work focus upon developing more sophisticated filters and learning algorithms, as these  Table 6: The runtime performance of individual bolts using random sift (randomly make sift decision with 50:50 split). Random sift was used to place greater processing load upon downstream bolts in the topology. Here process latency is the time taken to 'ack' a tuple after it is received. Note that 'acking' a tuple involves sending an acknowledgement to the transmitter of a tuple, so it knows it has been received. Execute latency is the time taken for the bolt to complete processing a tuple. Experiments run using c ratio = 0.1.
will likely greatly improve the pulsar recall rate. The algorithms installed at the bolts in the topology should also be studied in greater detail. There is scope to improve their runtime performance, and more crucially, filtering accuracy. We are currently working on a data generator that will assist us in such an investigation.