Single-pulse classifier for the LOFAR Tied-Array All-sky Survey
D. Michilli, J. W. T. Hessels, R. J. Lyon, C. M. Tan, C. Bassa, S. Cooper, V. I. Kondratiev, S. Sanidas, B. W. Stappers, J. van Leeuwen
MMNRAS , 1–10 (2017) Preprint 17 August 2018 Compiled using MNRAS L A TEX style file v3.0
Single-pulse classifier for theLOFAR Tied-Array All-sky Survey
D. Michilli, , (cid:63) J. W. T. Hessels, , R. J. Lyon, C. M. Tan, C. Bassa, S. Cooper, V. I. Kondratiev, , S. Sanidas, , B. W. Stappers, J. van Leeuwen , Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, 7990 AA, Dwingeloo, The Netherlands Astro Space Centre, Lebedev Physical Institute, Russian Academy of Sciences, Profsoyuznaya Str. 84/32, Moscow 117997, Russia Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Searches for millisecond-duration, dispersed single pulses have become a standard toolused during radio pulsar surveys in the last decade. They have enabled the discoveryof two new classes of sources: rotating radio transients and fast radio bursts. However,we are now in a regime where the sensitivity to single pulses in radio surveys is oftenlimited more by the strong background of radio frequency interference (RFI, whichcan greatly increase the false-positive rate) than by the sensitivity of the telescopeitself. To mitigate this problem, we introduce the Single-pulse Searcher (
SpS ). This isa new machine-learning classifier designed to identify astrophysical signals in a strongRFI environment, and optimized to process the large data volumes produced by thenew generation of aperture array telescopes. It has been specifically developed forthe LOFAR Tied-Array All-Sky Survey (LOTAAS), an ongoing survey for pulsarsand fast radio transients in the northern hemisphere. During its development,
SpS discovered 7 new pulsars and blindly identified ∼
80 known sources. The modular designof the software offers the possibility to easily adapt it to other studies with differentinstruments and characteristics. Indeed,
SpS has already been used in other projects,e.g. to identify pulses from the fast radio burst source FRB 121102. The softwaredevelopment is complete and
SpS is now being used to re-process all LOTAAS datacollected to date.
Key words: pulsars: general – surveys – methods: data analysis – methods: statistical
The first pulsar was discovered by recording its single pulsesat ∼ MHz using an aperture array (Hewish et al. 1968).In later studies, folding and Fourier-based techniques wereused to take advantage of the pulsar periodicity. Many pulsarobservations shifted to higher observing frequencies around . GHz, where the separation between pulsar signals andsky brightness is maximum for most of the pulsar popula-tion (Clifton & Lyne 1986). Furthermore, phased arrays weregenerally replaced by large single dishes, which remove thecomplexity of signal correlation permitting an increase intelescope sensitivity and bandwidth (Garrett 2012). How-ever, in recent years the increase in available computingpower makes it possible to build phased aperture array tele-scopes with sensitivities and bandwidths that outperform (cid:63)
E-mail: [email protected] traditional single dishes at low radio frequencies, offering alarger field-of-view (FoV) and more flexible instruments (vanHaarlem et al. 2013; Taylor et al. 2012; Tingay et al. 2013).This enables an exploration of a parameter space that iscomplementary to other searches: e.g., it is possible to detectsources having a spectrum steeper than the sky background(spectral index α ∼ − . , Mozdzen et al. 2017), which arelikely too faint to be detected at higher frequencies. In ad-dition, the larger FoV improves survey speed, and makesall-sky searches tractable.The most sensitive phased aperture array telescope todate is the LOw Frequency ARray (LOFAR, van Haarlemet al. 2013; Stappers et al. 2011). For example, its large col-lecting area already enabled detailed studies of archetypalsources, and the properties of the known pulsar population.The former is exemplified by the discovery of radio/X-raymode switching in PSR B0943+10 (Hermsen et al. 2013).Examples of the latter are the pulsar census results pre- © a r X i v : . [ a s t r o - ph . I M ] A ug D. Michilli et al. sented by Kondratiev et al. (2016) and Bilous et al. (2016).Beyond these known sources, many pulsars in our Galaxyremain undetected. It was recognised early on that LOFARhas great potential for discovering these (van Leeuwen &Stappers 2010). Pilot surveys placed limits on the occur-rence of fast transients at low frequencies, and discoveredthe first two radio pulsars with LOFAR (Coenen et al. 2014).We are now performing a full, sensitive survey of the north-ern hemisphere called the LOFAR Tied-Array All-Sky Sur-vey (LOTAAS, Coenen et al. 2014; Sanidas et al. 2018) .LOTAAS has already demonstrated its ability to find newpulsars using periodicity searches (with over 80 discoveriesto date) and in this paper we focus on discoveries madethrough single pulses. The long dwell time (1 hr) and largeFoV ( ∼ sq. deg. per pointing) of LOTAAS also make thesurvey potentially well placed to discover fast radio tran-sients, as long as they are not strongly affected by scatteringor dispersive smearing. An increasing issue for pulsar surveys is the presence ofradio-frequency interference (RFI) produced by several de-vices, which can mimic the behaviour of astrophysical sig-nals and limit survey sensitivities (e.g. Lyon et al. 2016).The large number of RFI detections makes it impracticalto visually inspect and follow-up all the detected signals.This is a worsening problem caused by the increasing num-ber of devices emitting radio waves and the improvements intelescope characteristics, such as sensitivity, dwell time andbandwidth. Therefore, this is a major challenge for next-generation radio telescopes and in particular the SquareKilometre Array (SKA, Ellingson 2004). An improvementis obtained by building telescopes in RFI-free zones, areaswhere the human presence is minimal. However, the radioemissions of air planes and satellites are still present.In order to lower the number of detections to be in-spected by eye, many automated classifiers have been devel-oped for pulsar surveys (see Lyon et al. 2016 and referencestherein for a summary). These automatic classifiers evolvedfrom simple heuristics and thresholds on S/N (e.g. Clifton &Lyne 1986) to semi-automated ranking algorithms (e.g. Leeet al. 2013). Also the graphical representation of the detectedsignals evolved, visualizing increasing information describingtheir parameters (e.g. Burgay et al. 2006). Machine learn-ing (ML) techniques began to be used to evaluate heuristicperformance, set threshold values and perform the classifi-cation (e.g. Eatough et al. 2010). Most recently, significantefforts are being spent in the selection of pulsar signals todeal with the large number of detections from new, sensitiveradio telescopes (e.g. Yao et al. 2016; Ford 2017; Bethapudi& Desai 2017).Here a branch of ML called statistical classification(Mitchell 1997) is used to filter to select astrophysical sig-nals. This requires first acquiring a large set of candidate ex-amples for which the ground truth origin or ‘class’ is known.When there are two classes under consideration, the classi-fication task is termed ‘binary’. In binary problems the tar-gets, i.e. astrophysical signals, belong to the positive class. The negative class describes all other examples (e.g. RFI ornoise). In either case the examples must be characterised viaone or more variables commonly known as ‘features’. Fea-tures are numerical or textual descriptors that summarise acandidate in some relevant way (e.g. S/N, pulse width, etc.).Features must be extracted by algorithms for each candi-date, and linked to their true class ‘labels’. When combined,this information forms what is known as a ‘training set’.Using supervised learning it is possible to ‘learn’ a mathe-matical function from this training set, that can automati-cally perform a similar mapping on new data. This processis known as ‘training’. The training process aims to splitthe training data into their respective classes, by using theinherent differences in the feature distributions to separatethem. The learning process is guided via a heuristic, mostoften error minimisation, that quantifies how many errorsare made. Each correct positive classification is known as aTrue Positive (TP), whilst an incorrect positive classifica-tion is known as a False Positive (FP). Similarly, negativeclassifications can be described in terms of True Negative(TN) and False Negative (FN) predictions. Together, thesefour outcomes form the so-called confusion matrix, used toassess how successfully an algorithm has learned the map-ping.For a good classification, it is essential to have featuresthat permit to separate the data into positive and negativeclasses. Therefore, specific algorithms must be designed toextract such features. Moreover, these algorithms need to befast and efficient in order to keep the computing time low.These algorithms can be designed by looking at the differ-ent properties of the members of the positive and negativeclasses. The classification success is usually determined dur-ing a ‘testing’ phase, conducted on an independent sampleof candidate examples. If the model learned during trainingperforms well during testing (produces few FPs and FNs) itcan be used to derive predictions for new previously unseendata. An algorithm will usually be successful if trained on alarge representative sample of data. Different heuristics canbe used to evaluate single features, such as the informationgain (also known as mutual information, Brown et al. 2012),a measure of the correlation between a feature and the tar-get variable (Lyon et al. 2016). Also, several metrics existto evaluate the performance of the whole set of features inclassifying the data. In this study, we made use of standardmetrics such as the false negative rate (FNR) and the falsepositive rate (FPR), which must be as low as possible for agood classification, the true positive rate (TPR, also knownas recall), the positive predictive value (PPV, also known asprecision), the accuracy (ACC), the G-Measure (G-M, thegeometric mean of recall and precision) and the F-Measure(F-M, the harmonic mean of recall and precision), whichmust be as high as possible for a good classification (e.g.Powers 2011). All these metrics assume values between zeroand one.
Soon after the initial discovery of pulsars, surveys began tak-ing advantage of the inherent periodicity of pulsar signals toimprove search sensitivity (Lorimer & Kramer 2004). Thisis usually achieved by performing a Fourier transform of thetime-series. However, this technique greatly decreases the
MNRAS000
MNRAS000 , 1–10 (2017) ingle-pulse Searcher sensitivity to sources whose emission is not regular over time(McLaughlin et al. 2006). Therefore, it has become standardprocedure to include single-pulse searches in pulsar surveysto avoid missing sources with large pulse amplitude vari-ations or a large null fraction (Lorimer & Kramer 2004).Moreover, two new classes of sources discovered in recentyears have created new interest in single-pulse searches: ro-tating radio transients (RRATs, McLaughlin et al. 2006) andfast radio bursts (FRBs, Lorimer et al. 2007). The formerclass is composed of pulsars whose emission is so sporadic intime that they are missed by periodicity searches. TypicalRRAT pulse rates range from as many as one per few sec-onds, up to one per several hours . The FRB class (Thorn-ton et al. 2013) is composed of extra-galactic radio flashes(Tendulkar et al. 2017). To date, only one has been observedto repeat (Spitler et al. 2016) and no periodicity has beendetected (Scholz et al. 2016).A typical search for single pulses is performed after de-dispersing the data collected by the telescope. This aimsto correct for the frequency-dependent delay induced by theinteraction of the radio waves from the source with free elec-trons present along the line of sight. The amount of disper-sive delay exhibited by a signal is proportional to the dis-persion measure (DM, which is the column density of freeelectrons). Since the DM of the source is unknown a priori ,it is necessary to de-disperse the signal at many trial values.The time-series resulting from the addition of all frequencychannels is then searched for single pulses. This is usuallyachieved by convolving each de-dispersed time-series witha top-hat function of variable width W . The properties ofthe function and the convolution are recorded every timethe resulting signal-to-noise ratio (S/N) is above a certainthreshold (see Lorimer & Kramer 2004 for a detailed discus-sion). The number of signal detections arising from RFI is par-ticularly large for LOTAAS because the LOFAR Core is ina region of high population density. In addition, the highsensitivity and long dwell time offered by the survey andthe large parameter space necessary to be searched at theselow frequencies increase the number of false detections. Forthese reasons, an automated classifier has been developed toclassify the periodic signals of LOTAAS (Lyon et al. 2016;Tan et al. 2018). It uses the advantages of statistical clas-sification to build a solid statistical framework for rejectingRFI.The presence of RFI is particularly problematic forsingle-pulse searches. In fact, as opposed to periodicitysearches, it is not possible here to filter out aperiodic signals.Usually, single-pulse classifiers take advantage of pulse shapein the frequency-time domain where, as opposed to RFI,astrophysical signals are expected to typically be broad-band and dispersed (Lorimer & Kramer 2004). As opposedto periodicity searches discussed earlier, only a few classi-fiers have been reported to specifically sift single-pulse de-tections. Keane et al. (2010) subtracted the non-dispersed http://astro.phys.wvu.edu/rratalog signal from their data (Eatough et al. 2009) and used spa-tial information from a multi-beam survey to discover tennew RRATs. Spitler et al. (2012) developed a multi-momenttechnique useful for quantifying the deviation present inthe pulse intensity at different frequencies. Karako-Argamanet al. (2015) presented RRATtrap , which uses manually-set thresholds to discriminate astrophysical signals basedon their S/N behaviour as a function of DM. Similarly,Deneva et al. (2016) developed
Clusterrank , which rejectsRFI instances that deviate from the theoretical relation be-tween signal strength, width and DM (Cordes & McLaughlin2003). Devine et al. (2016) studied different ML algorithmsapplied to their single-pulse classifier. For binary classifica-tion, a Random Forest (RF) trained on an oversampled set(RF over ) performed the best in their case.Here we present a new single-pulse classifier, Single-pulse Searcher ( SpS , Michilli & Hessels 2018), which isable to discriminate astrophysical signals from RFI withhigh speed and accuracy. In its current implementation(LOTAAS Single-pulse Searcher,
L-SpS ), it has been specif-ically designed to process LOTAAS data. LOTAAS obser-vations and data processing are presented in §
2; the
SpS classifier is presented in §
3; first discoveries are presented in §
4; and conclusions and future developments are discussedin § LOFAR is a radio telescope composed of thousands of an-tennas, which are grouped into stations that are distributedacross the Netherlands and other European countries (vanHaarlem et al. 2013; Stappers et al. 2011). LOTAAS ob-servations (Coenen et al. 2014) are performed using the sub-stations of the Superterp, a circular area with ∼ -m diameter where the station concentration is highest. TheHigh-Band Antennas (HBAs) are used to observe between − MHz. Signals from different stations are added bothcoherently (‘tied-array’ beams) and incoherently. The beamsare divided over three sub-array pointings (SAPs, which arepointing directions formed at station level). A total of 222 si-multaneous beams on sky are produced per pointing: inco-herent beams, a hexagonal grid of × coherent beams and × additional coherent beams that are scattered aroundthe central grid, and which can be pointed to known sourceswithin the FoV. The incoherent beams have a total FoV of ∼ deg , while the coherent beams have a total FoV of ∼ deg and a sensitivity ∼ √ = . times higher thanthat of the incoherent beams. Each LOTAAS pointing has a1-hr dwell time and a time resolution of . ms. The wholesurvey is divided into three passes and, upon survey com-pletion, each sky position will be covered three times by theincoherent beams and once by a coherent beam.Data are processed using a pipeline based on PRESTO (Ransom 2001) described in detail by Sanidas et al. (2018),which runs on the Dutch national supercomputer Cartesius .After using the rfifind algorithm to filter RFI in the time-frequency domain, the signal is incoherently de-dispersed https://userinfo.surfsara.nl/systems/cartesius MNRAS , 1–10 (2017)
D. Michilli et al. and frequency channels are added together using mpiprep-subband in order to form time-series data. Each time-seriesis searched for both periodic signals and single-pulse peaks.The latter is performed using single pulse search.py ,which convolves box-car functions of various lengths between . and ms. Single-pulse peaks with a S/N higher than5 are stored for further grouping and sifting.A total of ∼ DM trials is performed between DMvalues of – pc cm − with steps between DM trials of . – . pc cm − . Given a total of ∼ × time sam-ples per time-series, this implies a grid in the DM – timespace with × pixels, each one potentially containingan astrophysical signal, for each of the beams in everypointing. Considering only ideal random noise, this impliesan expected number of spurious detections ( ≥ σ ) for eachobservation of ∼ (Cordes & McLaughlin 2003). Instead, atypical observation produces ∼ - detections above σ ,demonstrating the huge impact that RFI and non-stationarynoise have on the single-pulse search. For comparison, thenumber of detections generated above σ is roughly an or-der of magnitude less and those above σ are normally halfof that. An automated algorithm capable of sifting thesedetections to identify the rare astrophysical signals is thusnecessary. L-SpS (Michilli & Hessels 2018) is a new sifting algorithmdesigned for single-pulse searches in a strong-RFI environ-ment, and specifically designed for the LOTAAS survey. Thesoftware uses ML techniques to differentiate RFI from as-trophysical signals, classify interesting signals and producediagnostic plots.
The aim of the program is to produce diagnostic plots for thevery best detections in a typical pulsar single-pulse search.Three steps are necessary to achieve this result.
The script single pulse search.py included in
PRESTO outputs information on each detected signal, specifically theDM, the time with respect to the beginning of the time-seriesand the selected width of the kernel function. These valuesare stored for each signal as one line in a text file. We defineone of such lines as an event . The
PRESTO -based pipelineis run for every LOTAAS pointing. For each beam,
L-SpS copies the events detected into the memory of the computerand stores them into a readily accessible HDF5 database .An impulsive signal having a large enough S/N, eitherRFI or astrophysical, will be detected in multiple time-seriesde-dispersed at nearby values. Therefore it will produce anumber of events close in time and DM. A broadband signalexhibits a decreasing S/N as the trial DM used moves fur-ther from the actual DM value because the pulse becomes increasingly smeared in time. This smearing is not symmet-ric and advances or delays the events for DM values higheror lower than the actual value, respectively. Therefore, theevents from a broadband signal will lie on a line in the DM –time space whose slope is given by ∆ t ∆ DM = k (cid:16) ν − m − ν − M (cid:17) , (1)where k = ( . ± . ) MHz pc − cm s is the disper-sion constant (Lorimer & Kramer 2004). This equation ishalf the value of the DM delay between the minimum ( ν m )and maximum ( ν M ) observing frequencies. In L-SpS , theslope defined by Eq. 1 is corrected so that events belongingto the same broadband signal are simultaneous.
All events close in time and DM are grouped together toform a pulse . The grouping is performed by a friends-of-friends algorithm (e.g. Huchra & Geller 1982; Press & Davis1982) designed to be highly efficient and able to process onemillion events in approximately half a second. Such an al-gorithm starts from the first event in the list and labels theclosest other event in the DM – time space to be part of thesame pulse if they are within a certain range. The thresh-olds on time and DM ranges between successive events wereset empirically to ms and DM trials (i.e. between . and pc cm − depending on the DM trial step size), respec-tively. The algorithm subsequently processes all the eventswith the same conditions. The characteristics of each pulse(i.e. time of arrival, DM, width, S/N) are derived from thebrightest event forming that pulse. Pulses formed by lessthan five events are considered spurious or too weak for sub-sequent analysis and therefore removed. Weak narrow-bandsignals, which produce pulses with constant S/N, are miti-gated by removing pulses with S/N < . . Also, pulses withDM < pc cm − are removed to avoid the contamination ofbroad-band RFI near DM = pc cm − . Tests are ongoing totry to lower this threshold on the DM. ML classification isthen applied to the pulses in order to discriminate RFI fromastrophysical signals, as discussed in § . pc cm − around the pulse. Also, these eventsmust have a S/N larger than half the pulse’s S/N in orderto only remove signals with a nearly constant strength overmany beams. If more than four additional beams containevents selected with those criteria, then the pulse is dis-carded. The spatial comparison is computationally expen-sive because of the many events to load and select. There-fore, it is implemented after the ML classifier has alreadyremoved a large fraction of the pulses. MNRAS000
All events close in time and DM are grouped together toform a pulse . The grouping is performed by a friends-of-friends algorithm (e.g. Huchra & Geller 1982; Press & Davis1982) designed to be highly efficient and able to process onemillion events in approximately half a second. Such an al-gorithm starts from the first event in the list and labels theclosest other event in the DM – time space to be part of thesame pulse if they are within a certain range. The thresh-olds on time and DM ranges between successive events wereset empirically to ms and DM trials (i.e. between . and pc cm − depending on the DM trial step size), respec-tively. The algorithm subsequently processes all the eventswith the same conditions. The characteristics of each pulse(i.e. time of arrival, DM, width, S/N) are derived from thebrightest event forming that pulse. Pulses formed by lessthan five events are considered spurious or too weak for sub-sequent analysis and therefore removed. Weak narrow-bandsignals, which produce pulses with constant S/N, are miti-gated by removing pulses with S/N < . . Also, pulses withDM < pc cm − are removed to avoid the contamination ofbroad-band RFI near DM = pc cm − . Tests are ongoing totry to lower this threshold on the DM. ML classification isthen applied to the pulses in order to discriminate RFI fromastrophysical signals, as discussed in § . pc cm − around the pulse. Also, these eventsmust have a S/N larger than half the pulse’s S/N in orderto only remove signals with a nearly constant strength overmany beams. If more than four additional beams containevents selected with those criteria, then the pulse is dis-carded. The spatial comparison is computationally expen-sive because of the many events to load and select. There-fore, it is implemented after the ML classifier has alreadyremoved a large fraction of the pulses. MNRAS000 , 1–10 (2017) ingle-pulse Searcher In every beam, positively classified pulses (occurring at dif-ferent times) are grouped into candidates according to theirproximity in DM. The maximum DM spread over which twopulses are considered to come from the same source is em-pirically set to . pc cm − , corresponding to 3–30 DM tri-als. Candidates in different beams within this DM range areconsidered to be produced by the same signal and only thebrightest is retained. The characteristics of each candidateare computed, i.e. beam number, average DM, cumulativeS/N, number of pulses detected and time of arrival if onlyone pulse is present, otherwise the period is calculated usingthe PRESTO routine rrat period . Candidates formed bya single pulse are considered only if they satisfy S/N > since weaker candidates would not be visible in the dynamicspectrum, meaning that their astrophysical nature can notbe verified. For the same reason, candidates formed by mul-tiple pulses must have a cumulative S/N > .A maximum of ten bright candidates (i.e. candidateswith the largest cumulative S/N) are processed per beam.Typically, a lower number of candidates per beam is pro-duced, as discussed in § ∼ less candidates. Diagnostic plots arethen generated for every selected candidate, as discussed in § Astrophysical signals of interest are different from RFI inthat they are usually broadband, dispersed and localized onthe sky, while normally the latter is either narrow-band ornot dispersed and often detected in multiple beams becausethe source is local. This implies that astrophysical signalsproduce pulses that peak in S/N in a certain beam and ata certain DM > pc cm − , while RFI will often have sim-ilar S/N in all the beams and constant S/N at differentDM trials if narrow-band, or peaked at DM = pc cm − if it is broadband. However, the real-world situation is oftenmore complicated because of statistical and non-stationarynoise (which tends to mask these differences especially forweak signals) as well as artefacts introduced by the tele-scope. Also, RFI can mimic a dispersive delay (e.g. Petroffet al. 2015), for example when multiple signals occur simul-taneously in the de-dispersed data. An astrophysical signalcan also be masked by simultaneous RFI. In addition, bothRFI and astrophysical signals can have complex structuresthat complicate their S/N curve as a function of DM. Fi-nally, some sources of interference like aeroplanes, radarsand signals reflected by the ionosphere can be beamed andthus appear localized on the sky.We make use of ML techniques in order to effectivelydifferentiate RFI and astrophysical signals and to have a sta-tistical foundation to evaluate the performance of our fea-tures. The statistical algorithm chosen to build the modeland perform the classification is the Gaussian-Hellinger VeryFast Decision Tree (Lyon et al. 2014), a tree learning algo- rithm (Mitchell 1997) based on the Very Fast Decision tree(VFDT; Hulten et al. 2001). The algorithm is designed todeal with imbalanced problems, where the target class, inour case astrophysical single-pulse events, are outnumberedby the instances we wish to reject (Lyon et al. 2016). Thedata mining tool WEKA was used to run this algorithmin order to evaluate classification performance, to build aclassification model and to perform the classification of newinstances. The other algorithms included in the standard WEKA distribution have been tested using the training setdescribed in the next subsection but none clearly outper-formed the VFDT.
A set of manually-labelled instances was selected to evalu-ate the features used and to build the classification model.A total of 35,063 instances of RFI were randomly chosenfrom various LOTAAS observations where no astrophysicalsources could be identified by eye. In addition, 18,003 pulseswere chosen from known pulsars selected in an equal num-ber of beams. All the pulses in these beams were includedin the training set if they were detected at the pulsar DMand were within . of the expected time of arrival giventhe pulsar period. The same conditions applied to randomDM and period values, in beams without any known pulsar,usually yielded no pulses to be selected. Given the costs as-sociated with labelling data, the possibility of human error,and because the ground truth labels can only be confirmedvia re-observing, it is possible for some training examples tobe incorrectly labelled. Though this is unlikely, a small pro-portion of training samples may be affected in this way. Suchmislabelled examples, whilst undesirable, have little impacton our results. This is because we used a sufficiently largecollection of labelled training samples describing all classes,which compensates. In any case, such mislabelled examplescan often help prevent a classifier from over-fitting, whichreduces real-world performance (Duda et al. 2000; Bishop2006). Using this selection process, also pulses affected bynoise or simultaneous RFI, which are desirable to be re-tained, were included. In this way, we aimed to reduce po-tential classification bias against the rare class. The distribu-tion of some parameters of RFI and astrophysical instancesin the training set is shown in Fig. 1. We developed different algorithms to extract features usedto select relevant signals. These algorithms were created byanalysing individual characteristics of RFI instances. At theend of the process, individual features had been evaluatedaccording to the value of their information gain. Featureswith a low information gain value were removed until a peakin the number of correctly classified instances was reached.Residual redundant features were subsequently removed viaan iterative method, manually deleting all the features oneby one and calculating the number of correctly classified in-stances for each configuration, until a peak was reached.Deneva et al. (2016) developed a classification algorithm MNRAS , 1–10 (2017)
D. Michilli et al.
Figure 1.
Comparison of the duration, S/N and DM of the RFI(black lines) and astrophysical (red lines) instances used in theLOTAAS training set. based on Eq. 12 of Cordes & McLaughlin (2003) that theyreported performing well. However, we could not use thesame algorithm to extract new features for
L-SpS becauseit was too slow in our implementation and required a largenumber of events per pulse.At the end of this process, the following five featureswere found to yield the highest number of correctly classifiedinstances. They are sorted in descending order according totheir information gain (Table 1):(i) Pulse width, W.(ii) Weighted mean of the pulse DM.(iii) Excess kurtosis of the width curve as a function ofDM. This curve is represented in red in Fig. 3h.(iv) Excess kurtosis of the S/N curve as a function of DM.This curve is represented in black in Fig. 3h.(v) Pulse S/N.Features (ii), (iii) and (iv) have been adapted from Tan et al.(2018).
The classifier’s ability to select astrophysical detections andto reject RFI instances is estimated using the manually-selected sample of pulses described in § . of instances in the LOTAAS set Table 1.
ML features selected for
L-SpS to separate astrophysicalpulses from RFI and their value of information gain. The subscripte indicates a single event within the pulse and σ is the standarddeviation. The rest of the symbols are defined in the text.Feature Inf. gain(i) W = W e ( max ( S/N e )) . (ii) DM = (cid:205) e DM e S/N e (cid:205) e DM e . (iii) k W = (cid:205) e ( DM e − DM ) W e σ ( W e ) (cid:205) e W e − . (iv) k S/N = (cid:205) e ( DM e − DM ) W e σ ( S/N e ) (cid:205) e S/N e − . (v) S/N = max ( S/N e ) . Table 2.
Confusion matrix for the ML classifier of
L-SpS appliedto the LOTAAS set. Predicted Predictedpulsar RFIActual pulsar TP = 17754 FN = 249Actual RFI FP = 355 TN = 34708 are correctly classified and the resulting confusion matrix isreported in Table 2.A manual inspection of the FP and FN instances re-turned by the classifier on the training set was performed.The vast majority of the mislabelled instances would havebeen misclassified even after visual inspection, with highprobability. It is possible that these instances were spuri-ous detections incorrectly labelled when the training set wasbuilt (see § As mentioned in §
2, the number of events that are storedat the end of the
PRESTO -based pipeline is on the order
MNRAS000
MNRAS000 , 1–10 (2017) ingle-pulse Searcher Table 3.
Metrics to evaluate the performance of different single-pulse classifiers. The ML algorithm of
L-SpS is applied to the LOTAAStraining set. The metrics used can assume any value between 0 and 1 and are described in the text. While the FNR and FPR assume lowervalues for a better classification, the rest of the metrics assume higher values for a better classification.
RRATtrap and
Clusterrank do not use ML and only rough estimations are available. The metrics are evaluated on different datasets and thus they are reported onlyfor reference (see main text for details). Values for RF over are from Devine et al. (2016), for RRATtrap are from Karako-Argaman et al.(2015) and for
Clusterrank are from Deneva et al. (2016).Classifier FNR FPR TPR ACC PPV G-M F-M
L-SpS .
014 0 .
010 0 .
986 0 .
989 0 .
980 0 .
983 0 . RF over .
282 0 .
011 0 . − − − . RRATtrap . . − − − − − Clusterrank . . − − − − − − )51015 S / N TP . . − )11121314 S / N FN . . − )6810 S / N FP . . . − )14 . . . S / N TN Figure 2.
S/N curves along DM of four pulses with the relativeclassification indicated on top. of - above σ for a typical LOTAAS observation. Byonly selecting one pulse for each group of events and by re-moving all the pulses formed by less than events, havingS/N < . or DM < pc cm − , as discussed in § ∼ . At this point, pulsesare selected using the ML classifier described in § ∼ pulses remain. The spatial comparison of differ-ent beams on the sky removes roughly half of these pulses.The grouping of pulses into candidates, the selection of thebrightest for each DM range and the thresholds on their cu-mulative S/N values discussed in § ∼ diagnostic files to inspect per observation. In total, L-SpS takes around 30 minutes to process one observation on oneCartesius node (powered by × cores . GHz Intel XeonE5-2690 v3 Haswell). Multimoment analysis could furtherreduce the number of RFI instances and tests are ongoing.
Despite the huge progress in ML classification in recentyears, a final human inspection is essential. Bright signalscan be easily identified by machines and by humans but forweak signals the classification becomes uncertain even afterhuman inspection. For this reason, we developed a set of di-agnostic plots to help visually identify astrophysical signalseven when they are buried in the noise. Once an observationhas been processed, a file containing the diagnostic plots isgenerated for each candidate. An example is shown in Fig. 3,which represents the discovery plot of PSR J0139+33. Forcomparison, its detection in a confirmation observation isshown in Fig. A1, while a diagnostic plot for an RFI candi-date is shown in Fig. A2.In the top part, four summary plots of detections in thebeam are present, along with meta data containing infor-mation on the beam where the candidate was detected. InFig. 3a, all the pulses positively classified in the beam areplotted with a size proportional to their S/N as a function oftime of arrival and DM. The pulses belonging to the candi-date shown are highlighted with stars and the ten brightestare numbered in order of decreasing S/N. Fig. 3b reports thenumber of positively classified pulses per DM range. Fig. 3cshows the distribution of S/N vs DM of the events formingthe pulses. Fig. 3d shows the cumulative S/N of the pulsesfor each DM range. The red lines indicate the candidate DM.In the bottom part, the sub-panels show informationrelated to the brightest pulse associated with the candidatesource. Meta data about the pulse is reported at the top-left. Fig. 3e is a close-up of the brightest pulse from Fig. 3a.It shows all the single events in the DM and time rangesplotted as black dots, including those considered to be noiseor RFI. Events forming the pulse are highlighted with redcircles, where the size is proportional to the S/N. The bluecross marks the most likely DM and time of the burst. Fig. 3fshows the time-series around the pulse time and DM. Redlines highlight the expected smearing of a broad-band sig-nal. Fig. 3g shows the profile of the pulse, i.e. the time-seriesaround the pulse’s time of arrival de-dispersed at the pulse’sDM. Fig. 3h represents the S/N (black line) and width (redline) of the events forming the pulse as functions of DM.Fig. 3i and j represent the pulse spectrum around the pulse’stime of arrival. These two plots are generated using water-faller.py . In Fig. 3i, the spectrum is de-dispersed at thepulse DM, down-sampled in frequency by a factor of 162 to https://github.com/plazar/pypulsar MNRAS , 1–10 (2017)
D. Michilli et al.
Figure 3.
Discovery plot of PSR J0139+33. Diagnostic plots like this are generated for each of the positively classified candidates fromthe LOTAAS survey. Top: summary plots of the pulses detected in the beam. Bottom: plots of the brightest pulse forming the candidate.The description of each sub-panel is given in §
16 sub-bands and smoothed in time with a boxcar functionof the same width as the pulse. The red triangle indicatesthe expected position of the signal. Fig. 3j shows the dis-persed data around the pulse, down-sampled in frequencyby a factor of 27 to 96 sub-bands and in time to 3 timesthe pulse width. The red curve represents the expected DM sweep of the signal. The RFI visible in Fig. 3i and j as hor-izontal stripes is the cause of the weak detections aroundthe pulse in Fig. 3e and of the short-duration events in thepulse tails in Fig. 3h. Fig. 3k reports the maximum S/N ofthe events detected in each of the coherent beams of theSAP, within the DM and time ranges indicated at the top.
MNRAS000
MNRAS000 , 1–10 (2017) ingle-pulse Searcher Table 4.
Pulsars discovered by early versions of
L-SpS in theLOTAAS survey. The level of pulse brightness variability duringthe observation is described qualitatively.Name Period (s) DM (pc cm − ) VariabilityPSR J0139+33 .
248 21 . ExtremePSR J0301+20 .
207 19 . LargePSR J0317+13 .
974 12 . LargePSR J0454+45 .
389 20 . SomePSR J1340+65 .
394 30 . SomePSR J1404+11 .
650 18 . SomePSR J1849+15 .
233 77 . Extreme
The S/N is normalized between 0 and 1 and represented asa hot-shades colour-scale. A white star indicates the beamwhere the candidate was detected and beams are numberedin blue. In this case, the signal is strongest in beam 73, wherethe pulse is detected. The complicated sensitivity patternand RFI is responsible for the (apparent) weaker detectionsin other beams.
The
L-SpS classifier has been developed and tested usingdata from the LOTAAS survey. The latest version of the pro-gram is presented in this paper and it will be used in a com-plete reanalysis of the survey data. Around 80 known pul-sars have been blindly identified by the different preliminaryversions of
L-SpS . Some of these were used to produce thetraining set described in § L-SpS — dozens have also been dis-covered using the periodicity searches and associated clas-sifier (Tan et al. 2018). A summary of their properties isreported in Table 4. PSR J0139+33 is a RRAT that was notdetectable in periodicity searches in the discovery pointing,nor in successive, targeted observations. PSRs J0301+20,J0317+13 and J1849+15 manifest a strong variability be-tween single pulses. For this reason, the first two were ini-tially detected through their bright single pulses and sub-sequently also found using a periodicity search in follow-upobservations. The remaining pulsars are bright enough tobe detected in both periodicity and single-pulse searches.The timing solutions of the sources presented here will bereported in future papers.
We have presented
L-SpS , a new classifier for searches of sin-gle radio pulses in the LOTAAS survey. It employs a ML al-gorithm to discriminate astrophysical signals from RFI, withhigh accuracy. During its development, the algorithm hasdiscovered 7 new pulsars and blindly identified ∼
80 knownsources. A full reprocessing of the LOTAAS data with thelatest version of
L-SpS , as presented here, is under way.Future improvements to the software include testing ofmultimoment analysis and development of additional fea-tures. Also, we only made use of binary classification, i.e.instances were divided into astrophysical and RFI. The use of multiclass classification (e.g. distinguishing broad-bandand narrow-band RFI) could improve the performance (Tanet al. 2018). In addition, a larger training set with positive in-stances better distributed (e.g. more high-DM pulsars) willbe used in future. Finally, deep learning techniques couldsignificantly improve the classification performance (Deng& Yu 2014). However, deep learning algorithms typically re-quire larger training sets. Therefore, they could possibly beused to reprocess the survey data when a sufficient numberof discoveries and re-detections is achieved.
Although
L-SpS has been designed specifically for theLOTAAS survey, efforts are ongoing to produce a more gen-eral software (
SpS , Michilli & Hessels 2018) that can bereadily adapted to other projects. The aim is to create aprogram that is user-friendly, simple to customize and ro-bust to different observation characteristics, in order to beeasily used in a generic study. This is achieved by designinga modular software where the different tasks discussed in § .At the time of writing, while the correct operation of SpS has been extensively tested, some of the features developedspecifically for the LOTAAS survey still need to be included,such as the capability to process multiple telescope beamsin parallel over different computer cores.Arguably the most critical part of the program is the fi-nal selection of astrophysical signals. In fact, due to the needfor a large data set of labelled detections, an ML analysis canbe difficult or impossible to perform in the case of small stud-ies. Therefore, a set of filters that do not rely on ML tech-niques was created for these situations. While such filtershave been designed to keep the false negative rate low, thefalse positive rate will be higher than in the ML approach,though still manageable for small projects. It is difficult toassess the general performance of these filters since they de-pend on the characteristics of the specific observations. Arough estimate is obtained from the study of FRB 121102with Arecibo. The features used typically reduce the num-ber of candidates by ∼ - %. Of the remaining candidates,typically ∼
25% were found to be real after visual inspection,though this fraction varied between 1% and 64% in differentobservations (depending on the severity of RFI).
ACKNOWLEDGEMENTS
This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.Computing time was provided by NWO Physical Sciences(project SH-242-15). The LOFAR facilities in the Nether-lands and other countries, under different ownership, are op-erated through the International LOFAR Telescope founda-tion (ILT) as an international observatory open to the globalastronomical community under a joint scientific policy. DM https://github.com/danielemichilli/SpS MNRAS , 1–10 (2017) D. Michilli et al. and JWTH acknowledge funding from the European Re-search Council under the European Union’s Seventh Frame-work Programme (FP/2007-2013) / ERC Starting Grantagreement nr. 337062 (“DRAGNET”). JWTH also acknowl-edges funding from an NWO Vidi fellowship. JvL acknowl-edges funding from the European Research Council un-der the European Union’s Seventh Framework Programme(FP/2007-2013) / ERC Grant Agreement n. 617199. Dif-ferent Python modules were used in this study, specifically numpy , matplotlib , pandas and astropy . REFERENCES
Bethapudi S., Desai S., 2017, preprint, ( arXiv:1704.04659 )Bilous A. V., et al., 2016, A&A, 591, A134Bishop C. M., 2006, Pattern Recognition and Machine Learning.SpringerBrown G., Pocock A., Zhao M.-J., Luj´an M., 2012, J. Mach.Learn. Res., 13, 27Burgay M., et al., 2006, MNRAS, 368, 283Clifton T. R., Lyne A. G., 1986, Nature, 320, 43Coenen T., et al., 2014, A&A, 570, A60Cordes J. M., McLaughlin M. A., 2003, ApJ, 596, 1142Deneva J. S., et al., 2016, ApJ, 821, 10Deng L., Yu D., 2014, Technical report, Deep Learning: Methodsand ApplicationsDevine T. R., Goseva-Popstojanova K., McLaughlin M., 2016,MNRAS, 459, 1519Duda R. O., Hart P. E., Stork D. G., 2000Eatough R. P., Keane E. F., Lyne A. G., 2009, MNRAS, 395, 410Eatough R. P., Molkenthin N., Kramer M., Noutsos A., KeithM. J., Stappers B. W., Lyne A. G., 2010, MNRAS, 407, 2443Ellingson S. W., 2004, Experimental Astronomy, 17, 261Ford J. M., 2017, PhD thesis, Nova Southeastern UniversityGarrett M. A., 2012, in From Antikythera to the SquareKilometre Array: Lessons from the Ancients. p. 41( arXiv:1211.6455 )Hermsen W., et al., 2013, Science, 339, 436Hewish A., Bell S. J., Pilkington J. D. H., Scott P. F., CollinsR. A., 1968, Nature, 217, 709Huchra J. P., Geller M. J., 1982, ApJ, 257, 423Hulten G., Spencer L., Domingos P., 2001, in Proceedings of theSeventh ACM SIGKDD International Conference on Knowl-edge Discovery and Data Mining. KDD ’01. ACM, New York,NY, USA, pp 97–106, doi:10.1145/502512.502529, http://doi.acm.org/10.1145/502512.502529
Karako-Argaman C., et al., 2015, ApJ, 809, 67Keane E. F., Ludovici D. A., Eatough R. P., Kramer M., LyneA. G., McLaughlin M. A., Stappers B. W., 2010, MNRAS,401, 1057Kondratiev V. I., et al., 2016, A&A, 585, A128Lee K. J., et al., 2013, MNRAS, 433, 688Lorimer D. R., Kramer M., 2004, Handbook of Pulsar AstronomyLorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J.,Crawford F., 2007, Science, 318, 777Lyon R., Brooke J., Knowles J., Stappers B., 2014, 2014 22ndInternational Conference on Pattern Recognition (ICPR), 00,1969Lyon R. J., Stappers B. W., Cooper S., Brooke J. M., KnowlesJ. D., 2016, MNRAS, 459McLaughlin M. A., et al., 2006, Nature, 439, 817Michilli D., Hessels J. W. T., 2018, SpS: Single-pulse Searcher,Astrophysics Source Code Library (ascl:1806.013)Mitchell T. M., 1997, Machine Learning, 1 edn. McGraw-Hill,Inc., New York, NY, USA Mozdzen T. J., Bowman J. D., Monsalve R. A., Rogers A. E. E.,2017, MNRAS, 464, 4995Petroff E., et al., 2015, MNRAS, 451, 3933Powers D. M. W., 2011, Journal of Machine Learning Technolo-gies, 2, 37Press W. H., Davis M., 1982, ApJ, 259, 449Ransom S. M., 2001, PhD thesis, Harvard UniversitySanidas S., et al., 2018, in prep.Scholz P., et al., 2016, ApJ, 833, 177Spitler L. G., Cordes J. M., Chatterjee S., Stone J., 2012, ApJ,748, 73Spitler L. G., et al., 2016, Nature, 531, 202Stappers B. W., et al., 2011, A&A, 530, A80Tan C. M., Lyon R. J., Stappers B. W., Cooper S., HesselsJ. W. T., Kondratiev V. I., Michilli D., Sanidas S., 2018, MN-RAS, 474, 4571Taylor G. B., et al., 2012, Journal of Astronomical Instrumenta-tion, 1, 1250004Tendulkar S. P., et al., 2017, ApJ, 834, L7Thornton D., et al., 2013, Science, 341, 53Tingay S. J., et al., 2013, Publ. Astron. Soc. Australia, 30, e007Yao Y., Xin X., Guo P., 2016, in 2016 12th International Con-ference on Computational Intelligence and Security (CIS). pp120–124, doi:10.1109/CIS.2016.0036van Haarlem M. P., et al., 2013, A&A, 556, A2van Leeuwen J., Stappers B. W., 2010, A&A, 509, A7
APPENDIX A: RFI DIAGNOSTIC PLOT
Examples of two diagnostic plots generated by
L-SpS arepresented. Fig. A1 reports the detection of PSR J0139+33during its confirmation observation. Fig. A2 shows an exam-ple of a typical diagnostic plot for a candidate produced byRFI. The explanation of the sub-plots is given in § This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000
L-SpS arepresented. Fig. A1 reports the detection of PSR J0139+33during its confirmation observation. Fig. A2 shows an exam-ple of a typical diagnostic plot for a candidate produced byRFI. The explanation of the sub-plots is given in § This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000 , 1–10 (2017) ingle-pulse Searcher Figure A1.
Example of diagnostic plots of the confirmation observation of PSR J0139+33. The same plots as Fig. 3 are represented andthey are discussed in the text. The number of beams visible in panel k is different for confirmation observations. The strong interferencevisible as horizontal stripes in panels i and j did not prevent to detect the bright pulses from the RRAT.MNRAS , 1–10 (2017) D. Michilli et al.
Figure A2.
Example of diagnostic plots for an RFI candidate. The same plots as in Fig. 3 are represented and they are discussed in thetext. MNRAS000