Deep learning approach for identification of HII regions during reionization in 21-cm observations
Michele Bianco, Sambit. K. Giri, Ilian T. Iliev, Garrelt Mellema
MMNRAS , 1–15 (2021) Preprint 16 February 2021 Compiled using MNRAS L A TEX style file v3.0
Deep learning approach for identification of
H ii regions duringreionization in 21-cm observations
Michele Bianco, , ★ Sambit. K. Giri, , Ilian T. Iliev , Garrelt Mellema Astronomy Centre, Department of Physics & Astronomy, Pevensey III Building, University of Sussex, Falmer, Brighton, BN1 9QH, United Kingdom The Oskar Klein Centre, Department of Astronomy, StockholmUniversity, AlbaNova, SE-10691 Stockholm, Sweden Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The upcoming Square Kilometre Array (SKA-Low) will map the distribution of neutral hydrogen during reionization, andproduce a tremendous amount of 3D tomographic data. These images cubes will be subject to instrumental limitations, such asnoise and limited resolution. Here we present
SegU-Net , a stable and reliable method for identification of neutral and ionizedregions in these images.
SegU-Net is a U-Net architecture based convolutional neural network (CNN) for image segmentation. Itis capable of segmenting our image data into meaningful features (ionized and neutral regions) with greater accuracy comparedto previous methods. We can estimate the true ionization history from our mock observation of SKA with an observation time of1000 h with more than 87 per cent accuracy. We also show that
SegU-Net can be used to recover various topological summarystatistics, such as size distributions and Betti numbers, with a relative difference of only a few per cent. These summary statisticscharacterise the non-Gaussian nature of the reionization process.
Key words: cosmology: radio – early universe, reionization – deep learning, image processing
The Epoch of Reionization (EoR) is a period of great importancein the study of structure formation and evolution in the Universe.During this period the predominately cold and neutral intergalac-tic medium (IGM) transitioned to a hot and ionized state due tothe appearance of the first luminous cosmic sources. These sources,which may have been star forming galaxies and quasi-stellar objects(QSOs), produced the ionizing photons which over a period of ap-proximately 500 million years completed the re-ionization of theUniverse. (Furlanetto et al. 2006; Zaroubi 2012; Ferrara & Pandolfi2014).This period is one of the least understood epochs in the historyof the Universe, due to the lack of direct observations. Indirect con-straints have been put on the reionization process with observationof Lyman- 𝛼 forest (e.g. Fan et al. 2006; McGreer et al. 2011, 2014),number density of Lyman- 𝛼 emitters (e.g. Ota et al. 2008; Ouchi et al.2010; Robertson et al. 2015), high-z quasar spectra (e.g. Schroederet al. 2013; Totani et al. 2016; Davies et al. 2018; Greig et al. 2019)and Thomson scattering optical depth of cosmic microwave back-ground (CMB) photons (e.g. Komatsu et al. 2011; Planck Collabo-ration et al. 2016).The ground state of neutral hydrogen atom can produce a signalthrough spin-flip transition, which is known as the 21-cm signal. Thissignal will be a unique signature of EoR (e.g. Madau et al. 1997;Furlanetto et al. 2006). When observed, the 21-cm signal would haveredshifted to radio band of the electromagnetic spectrum. Various ★ Contact e-mail: [email protected] radio experiments, such as Low Frequency Array (LOFAR; e.g.van Haarlem et al. 2013), Murchison Widefield Array (MWA; e.g.Tingay et al. 2013) and the Hydrogen Epoch of Reionization Array (HERA; e.g. DeBoer et al. 2017), have been trying to detect thissignal. Recently these facilities have provided useful upper limits onthe 21-cm power spectrum (e.g. Mertens et al. 2020; Trott et al. 2020)that have been used to learn the properties of reionization (e.g. Gharaet al. 2020; Mondal et al. 2020; Greig et al. 2020a,b).The 21-cm signal during EoR will be highly non-Gaussian andtherefore the power spectrum will not give the full characterisa-tion (e.g. Mellema et al. 2006a; Ichikawa et al. 2010; Watkinson& Pritchard 2015; Majumdar et al. 2018; Giri et al. 2019b). In thecoming years, the Square Kilometre Array (SKA) will be build.The low frequency component of the SKA will be sensitive enoughto detect the 21-cm signal produced during EoR and create imagesof its distribution on the sky (Mellema et al. 2013; Wyithe et al.2015; Koopmans et al. 2015). These images contain more informa-tion about our Universe as the detection of the signal at differentobserved frequencies depict the distribution of neutral hydrogen at agiven time during the EoR.SKA-Low will observe a sequence of such 21-cm images fromdifferent redshifts will constitute a three-dimensional set of data,which is known as tomographic dataset, and the evolution of the21-cm signal can be seen along the redshift axis. See for example http://reionization.org/ https://skatelescope.org © a r X i v : . [ a s t r o - ph . I M ] F e b M. Bianco et al.
Giri (2019) for more description about tomographic 21-cm images.Reionization process is driven by growing bubbles of H ii regions(e.g. Furlanetto et al. 2004). As the sources of ionizing photons willreside inside these bubbles, observing these bubbles and its evolutionwill be interesting. Numerous studies have provided various methodsto detect and study properties of H ii bubbles (e.g. Datta et al. 2007;Zackrisson et al. 2020; Mason & Gronke 2020). We can also studythe properties of reionization with 21-cm images (Giri et al. 2018a,2019a). However, tomographic images from SKA-Low will be proneto instrumental limitations, such as noise, limited resolution andforeground contamination (e.g. Koopmans et al. 2015; Ghara et al.2016). In the field of image processing, the methods that can classifyobjects or features in images into meaningful segments are knownas ‘image segmentation’. Giri et al. (2018b) implemented an imagesegmentation method to classify neutral and ionized regions in 21-cmimages in the presence of instrumental limitations and demonstratedthat we can derive properties of reionization from such observation.Artificial intelligence (AI) and deep learning methods are capableof learning patterns in image data and identifying interesting regions.Image segmentation is a common problem in the field of imageprocessing, aiming to learn visual patterns from input image andidentify objects classes that make up that image. This technology isquite popular in the field of data analysis and has been applied tostudy objects with complex visual form contained in big data (Longet al. 2014). In recent years, several papers made use of machinelearning techniques also in different problems of astrophysics (e.g.Lee 2019; Giri et al. 2020b; Yoshiura et al. 2020; Chen et al. 2020) andcosmology (e.g. Jeffrey et al. 2020; Sadr & Farsian 2020; Guzman &Meyers 2021). In the case of reionization, several of these methodsare aimed to either remove foreground emission (Li et al. 2019;Makinen et al. 2020; Villanueva-Domingo & Villaescusa-Navarro2020), emulate reionization simulations (e.g. Kern et al. 2017; Schmit& Pritchard 2018; Jennings et al. 2018; Cohen et al. 2019; Ghara et al.2020) or constrain reionization history (e.g. Chardin et al. 2019;Mangena et al. 2020) and its astrophysical inputs (e.g. Sullivan et al.2017; Gillet et al. 2019; Hassan et al. 2020).In this work, we present a new approach for the identification ofthe distribution of H ii regions in 21-cm images using a deep learningmethod named U-shaped convolutional neural network (U-Net), thatis especially designed for image segmentation and feature extraction(Ronneberger et al. 2015). In our case, we adapt this network forprocessing our image data, which are mock observations of the 21-cm signal during the EoR. The method will segment the images intoionized and neutral regions. We name this framework
SegU-Net .This paper is organized as follows. In § 2 we present how wegenerate the simulated data set used for this work. In § 3 we describethe design of our neural network, including the error estimation.In § 4 we discuss its application to our simulated SKA-Low datasets, considering a range of summary statistics such as the meanionization fraction, power spectra and topological quantities such assize distributions and Betti numbers. In § 5 we test our frameworkon various instrumental noise levels and in § 6 we test it on a data setproduced from a fully numerical reionization simulation. We discussand summarize our conclusions in § 7.
For any deep learning based method, we need a data set containingsample of the possible scenarios, known as the training set. In § 2.1,we describe the reionization simulation code that we use to create thetraining set. The observable for radio telescopes observing the 21-cm
Table 1.
The parameters used in this study to model the telescope properties.Parameters ValuesSystem temperature 60 ( 𝜈 ) − . KEffective collecting area 962 m Declination -30 ◦ observation hour per day 6 hoursSignal integration time 10 seconds signal is defined in § 2.2. Finally in § 2.3 we give the methodologywe use to mimic the observations expected with SKA-Low. In order to train our network, we require a large set of simulationsthat represent the 21-cm radio signal, for a wide range of redshiftduring reionization and different astrophysical properties. To do sowe employ py21cmFAST , the
Python wrapped version of the semi-numerical cosmological simulation code (Mesinger et al.2011; Murray et al. 2020). The code computes the evolution of thematter density field using the Zel’dovich approximation (Zel’Dovich1970). The ionization field and the corresponding 21-cm differentialbrightness temperature are then calculated from the matter densitydistribution based on the excursion set formalism (Furlanetto et al.2004; Mesinger & Furlanetto 2007), which consider a region to beionized when the fraction of collapsed matter fluctuation exceed amass threshold. The ionization fraction 𝑥 HII ( 𝑟𝑟𝑟 ) at a position 𝑟𝑟𝑟 isgiven as, 𝑥 HII ( 𝑟𝑟𝑟 ) = (cid:40) 𝑓 coll ≥ / 𝜁 ion 𝜁 𝑖𝑜𝑛 is the ionizing efficiency of high redshift galaxies and 𝑓 coll ( 𝑅 s , 𝑀 min ) is the fraction of collapsed matter within radius 𝑅 s that can form haloes with mass greater than 𝑀 min . 𝑓 coll is calculatedat every pixel varying 𝑅 s within 0 and 𝑅 mfp . The maximum valueof 𝑓 coll is used in Equation 1. 𝑅 mfp is the mean free path length forionizing photons in the IGM.The cosmological parameters considered in this work are based onWMAP 5 years data observation (Komatsu et al. 2009) and consistentwith Planck Collaboration (2019) results. We assume a flat Λ CDMcosmology with the following parameters, Ω Λ = . Ω 𝑚 = . Ω 𝑏 = . 𝐻 = 𝑘𝑚 𝑠 − 𝑀 𝑝𝑐 − , 𝜎 = . 𝑛 𝑠 = . Radio interferometry based telescopes record the differential bright-ness temperature 𝛿𝑇 b while observing the redshifted 21-cm signal. 𝛿𝑇 b depends on position on the sky 𝑟𝑟𝑟 and redshift 𝑧 and can be givenas (e.g. Mellema et al. 2013), 𝛿𝑇 b ( 𝑟𝑟𝑟, 𝑧 ) ≈ 𝑥 HI ( 𝑥𝑥𝑥, 𝑧 ) (cid:0) + 𝛿 b ( 𝑟𝑟𝑟, 𝑧 ) (cid:1) (cid:18) + 𝑧 (cid:19) (cid:18) − 𝑇 CMB ( 𝑧 ) 𝑇 s ( 𝑟𝑟𝑟, 𝑧 ) (cid:19)(cid:18) Ω b . ℎ . (cid:19) (cid:18) Ω m . (cid:19) − mK , (2)where 𝑥 HI , 𝛿 b , 𝑇 CMB and 𝑇 s are neutral fraction, baryon densitycontrast, CMB temperature and spin temperature respectively.Previous studies have shown that our Universe will be heatedbefore reionization begins (e.g. Pritchard & Furlanetto 2007; Rosset al. 2017, 2019). Therefore we assume 𝑇 s (cid:29) 𝑇 CMB throughout this
MNRAS , 1–15 (2021) Figure 1.
Top left, the hydrogen neutral fraction at simulation resolution. Green lines contours indicated the separation between neutral and ionized region,smoothed with a baseline of 𝐵 = 𝑘𝑚 . Bottom right, an example of a smoothed box slice used during the network training. To the differential brightnesstemperature at simulation resolution (Bottom left), we subtracted the mean signal in the frequency axis. We then added the Gaussian noise for 1000 hoursobserved time (Top right) and smoothed the cube with the same baseline. Black solid line shows the same contour as in top left panel. work, which is known as the spin saturated approximation and isrelevant at low redshift 𝑧 (cid:46)
12. In spin saturated approximation sce-nario, the differential brightness signal will be observed in emission( 𝛿𝑇 b ≥ ) . Thus the positions on the sky with 𝛿𝑇 b = In order to train
SegU-Net for application to actual observations,we need a training set of mock observations. We create these mockobservations by simulating the 𝛿𝑇 b using the methods described inprevious sections and adding instrumental effects, such as absence ofzero baseline, limited resolution and noise. We follow the methods inGhara et al. (2016) and Giri et al. (2018b) for mimicking the expectedeffects of SKA1-Low.We consider a simulation volume of (
256 Mpc ) and an intrinsicresolution of Δ 𝑥 = Δ 𝜃 = .
777 arcminand a frequency depth of Δ 𝜈 = .
124 MHz along the line of sightat 𝑧 =
7. As an example, in Figure 1, we show a coeval cube sliceof the neutral fraction field and 𝛿𝑇 b field in top left and bottom leftpanels respectively. These slices are taken from the epoch when theuniverse was about 50 per cent ionized. For each 𝛿𝑇 b coeval cube,we assume one axis as the line of sight or frequency direction and subtract the mean signal from each frequency channel, such thatthis could be considered as a sub-volume from the 3D tomographicdata set. We consider this simulation as our reference throughoutthe results analysis in §4, its astrophysical parameters are shown inTable 2.We simulate the instrumental noise using the method given inGiri et al. (2018b) and is implemented in Tools21cm (Giri et al.2020a). We provide our assumption for telescope setup in Table 1.In top right panel of Figure 1, we show a slice of from the simulatednoise cube produced with 1000 hours observation with SKA1-Lowat simulation resolution. When we add this noise to our simulatedsignal at the simulation resolution, we cannot see any feature of thesignal as the noise is several orders of magnitude higher than thesignal. Therefore we degrade the signal with noise to a resolutioncorresponding to a maximum baseline 𝐵 = 𝑥 th = . A python package for EoR simulations analysis. https://github.com/sambit-giri/tools21cm
MNRAS000
MNRAS000 , 1–15 (2021)
M. Bianco et al.
Figure 2.
Schematic representation of
SegU-Net network architecture. Orange arrow indicates 2D convolutional layer, followed by batch normalization andReLU activation. Pooling operation followed by dropout layer are in green arrow. Blue arrow indicates up-sampling layer by transposed 2D convolutional layerand with a red arrow the closing layer, a 2D convolution followed by a sigmoid activation function. The descending path on the left side, divide the resolutionof the image, after each pooling operation, and double the channel dimension after each convolution. On the other hand, the expansion path double the spatialdimension at each up-sampling operation and decrease the channel dimension after concatenation with its counterpart layer, in the descending path.
Table 2.
Astrophysical parameters used for our fiducial simulation.Parameters Values 𝜁 . 𝑅 mfp .
861 Mpc 𝑇 minvir . × K neutral/ionized regions. Then we over-plot the boundaries of theseionized regions, the neutral fraction slice and signal slice in top-leftand bottom-right panels of Figure 1 respectively. For our training set, we randomly sampled the astrophysical sim-ulation parameters by a normal distribution, such that the ioniz-ing efficiency of high-redshift galaxies 𝜁 is sampled with N ∼( . , ) , the mean free path of ionizing photons 𝑅 mfp with N ∼ ( . , ) and the (logarithmically-spaced) minimumvirial temperature for halos to host star-forming galaxies 𝑇 minvir with N ∼ ( . , . ) . The redshift is randomly sampled with an uni-form distribution U ∼ [ , ] . The initial conditions of the cosmo-logical density field are changed for each simulation. This helps usavoid the impact of cosmic variance on our trained model. With thelist of all the parameter values, we produced 10,000 mock observa-tions of 21-cm signal. Out of these mock observations, we use 15per cent as the so-called network validation set. This validation set isused during the training method to provides an unbiased evaluationof the network model fit.Eventually we will use SegU-Net on actual 21-cm image obser-vations. Here we rely on an additional 300 mock observations as thetesting or prediction set. Just as for the training set, the parametervalues are randomly chosen. We call this the ‘random’ testing set.The training process is blind to the prediction set.Apart from the above testing set, we create an additional simulationwith fixed values of astrophysical parameters (given in Table 2). Wehave chosen these values such that between 𝑧 = 𝑥 VHI ≈ . SegU-Net ’s capability to capture the evolutionof structures and recover the binary field from untrained data in § 4.
Here we describe our machine learning method for identifying ion-ized and neutral regions in noisy 21-cm images and its predictionuncertainty estimation in § 3.1 and § 3.2 respectively.
SegU-Net
Our segmentation network is based on the U-Net framework first in-troduced by Ronneberger et al. (2015). U-Net consists of two likewisesymmetric paths, an encoder operator that contracts the image anda decoder operator that expands the extracted features. The encodercorresponds to a classical convolutional neural network (CNN). Thepurpose of this CNN is to reduce the size of the input image in such away that only information of the most interesting features remains. Aseries of concatenated convolution operations (layers) return a low di-mensional latent space (or latent vector) which contains informationabout these extracted features. We show a schematic representationof the U-Net in Figure 2. The left part of the U-shape in the diagramand the bottom layer represent the encoder and the low dimensionallatent space respectively. For a detailed discussion of CNNs, we referthe reader to Mehta et al. (2019), and for examples of employingCNNs to infer cosmological and astrophysical parameters in the con-text of reionization to Gillet et al. (2019) and Hassan et al. (2020). Inour case, the information contained in the latent space (or latent vec-tor) of U-Net (bottom layer) is expanded by a decoder into a binarymap of the same size as the input image. The right part of the U-shape of the diagram in Figure 2 represents the decoder. The decodergradually increases the spatial resolution of the latent vector with anup-sampling operation (transposed convolution), until we obtain the https://github.com/micbia/SegU-Net MNRAS , 1–15 (2021) Figure 3.
Slice comparison of the binary field (in blue ionized region, and in red neutral) recovered by the Super-Pixel method (left) and our neural network(center). On the right panel, the per-pixel uncertainty calculated in
SegU-Net . Green lines indicates the true separation between ionized/neutral regions, derivedby the simulation neutral hydrogen distribution. same dimension of the input image. After each up-sampling step, theoutput is combined with the corresponding encoder layer with thesame dimension. We illustrate this further in Appendix B with anexample.Even though each of our image data is 3D,
SegU-Net is trainedon 2D slices. We identify structures in 3D image data by running onevery slice along the third axis. We tested and found that the methodis not sensitive to the choice of the third axis. When compared toa neural network trained on 3D data, we found that our approach iscomputationally less expensive without loss of accuracy. Thereforethe U-Net architecture described in this work will be for 2D imagedata.The structure of the encoder layers consists of two convolutionalblocks followed by a 2D max-pooling layer (
MaxPool ) of size 2x2and a 5 per cent rate dropout layer (
Drop ), a regularization tech-nique that randomly shut down a portion of the layer neurons toavoid over-fitting (Hinton et al. 2012; Srivastava et al. 2014). Theconvolutional block (
ConvBlock ) consists of a 2D convolution layer(
Conv2D ) with 3x3 kernel size. We add a layer that normalizes theprevious input layer over the batch sample to avoid over-fitting ( BN )(Ioffe & Szegedy 2015) and as an activation function we employa Rectified Linear Unit ( ReLU ) activator (Jarrett et al. 2009; Glo-rot et al. 2011),
ConvBlock=Conv2D+BN+ReLU . This layer structureis repeated for a total of four levels (
Encoder-Level ) and at eachstep the dimension of the input image is halved by a pooling oper-ation. The number of feature channels are doubled by the convolu-tional layer,
Encoder-Level=2*ConvBlock+MaxPool+Drop . Thedecoder structure is somewhat similar to the encoder, we replacethe pooling operation with a transposed 2D convolution (
TConv2D )(Dumoulin & Visin 2016; Zeiler & Fergus 2013), that has anopposite scaling effect on the resolution and channel size. Thislayer output is then concatenated ( CC ) with the corresponding en-coder level, in order to preserve the features extracted in the con-tracting path. This step is followed by a dropout layer, Decoder-Level=TConv2D+CC+Drop+2*ConvBlock . The final output con-sists of a 2D convolutional layer followed by a sigmoid activation.Our network has a total of 23 2D convolutional layers distributed onfour down- and up-sampling scaling levels and a bottom layer, for atotal of approximately 2.5 million trainable parameters. In Figure 2,we show our best performing network and label the shape of output from each intermediate hidden layer of this network. In Appendix Aand B we discuss this further.During our training process, the hyperparameters of the networkare learnt by minimizing a loss function. We employ the balancedcross-entropy (BCE) (Salehi et al. 2017), L( y , ˆy ) = − N ∑︁ i = ( 𝛽 y i log ( ˆy i ) + ( − 𝛽 )( − y i ) log ( − ˆy i )) (3)where 𝑦 𝑖 ∈ {
0; 1 } is the pixel-wise ground truth, ˆ 𝑦 𝑖 the predictedvalue, 𝑁 the batch size, which is our case is of size 32 and the param-eter 𝛽 = (cid:205) 𝑁𝑖 = 𝑦 𝑖 is the average volume neutral fraction of the batch.In our context, at early/late stage of reionization the statistical weightof the ionized/neutral pixels are under-represented. This situation isknown in data science as a problem affected by “class unbalanced”data. To deal with this we use the above loss function which hasbeen shown to be well suited for segmentation problems that areaffected by class unbalanced data (Cui et al. 2019). We further usedthe Adaptive Moment Estimator Adam (Kingma & Ba 2014), an op-timized stochastic gradient descent algorithm for error minimization.The initial learning rate, the step size of the rate of convergence thatminimizes the loss function, is set to 10 − . We run the network on 2GPUs, and it took approximately 1,500 core-hours. SegU-Net
One of the main drawbacks of machine learning is that it is unableto quantify uncertainties and confidence intervals for its predictions,and only recently attempt have been made (Charnock et al. 2020;Hortúa et al. 2020) to include error estimation. However, this has notyet been generally implemented for U-Nets. Additionally, if not welloptimized, neural networks are prone to over-fitting and tend to bebiased. We therefore have developed an error estimation procedureto be used during the prediction process. This gives our networkadditional power and provides a pixel by pixel error map.Image manipulations, such as zooming, shifting along an axis,flipping axes and rotation along an axis, are commonly performedon 2D or 3D image training data to increase the number of samples(Simonyan & Zisserman 2015; Szegedy et al. 2015). This techniqueis known as time-test augmentation (TTA) of data (Perez & Wang
MNRAS000
MNRAS000 , 1–15 (2021)
M. Bianco et al.
Figure 4.
On the left, the Matthews correlation coefficient 𝑟 𝜙 of the recovered binary field for the prediction set, against its volume averaged neutral fraction.Error-bar indicates the network confidence interval and colors indicates the redshift of the simulated coeval cube. On the inset panel, we show the distributionof the training set (blue step-line) against the volume average neutral fraction. On the right, we compare the same correlation coefficient for recovers performedfor the fiducial simulation with our neural network (blue circle line) and Super-Pixel method (orange square line). We also include the result on the test on C RAY simulation, from left to right, redshift 𝑧 = . , . , . , .
40 corresponding to a volume averaged neutral fraction of 𝑥 VHI = . , . , . , . ∼ SegU-Net . Each of the recovered binary fields is transformed back.We calculate the final result as the average of these fields and theper-pixel standard deviation as an estimate of the error for each pixel.An example of the pixel per pixel error map can be seen in Figure 3(right-most panel). We will discuss this figure further in § 4.1. Thissimple method provides our neural network with an error estimationfor each labeled pixel.
Once the network is trained we want to estimate how well it recoversthe binary field from noisy 21-cm images. To do so we include in ouranalysis the up-to-date Super-Pixel method presented in Giri et al.(2018b). In previous studies, this method has been shown to be quiteefficient in extracting useful summary statistics from noisy 21-cmimages (Giri et al. 2019a; Giri & Mellema 2020).
To start, we show a visual comparison of slices in Figure 3. We com-pare the predicted binary field recovered by the Super-Pixel method(left-most panel) and
SegU-Net (central panel) with the ground truth(green contours in both panels). The ground truth is boundary of ion-ized regions which are extracted from the simulation neutral fractionfield at same resolution by putting a threshold of 0.5. The red and bluepixels represent neutral and ionized pixels respectively. In the right-most panel we show the pixel-error estimated from
SegU-Net witha color-bar. The error is determined by calculating the standard-deviation of the same pixel from the different version of the samemock observation produced with TTA (see § 3.2).
Table 3.
Summary of the Matthews correlation coefficient score (in per cent)of our two test sets for the two feature identification methods.
SegU-Net
Super-Pixelredshift random set fiducial random set fiducial 𝑧 ≤ .
75 88.9% 91.7% 63.7% 62.6% 𝑧 ≥ .
25 85.3% 90.1% 60.7% 71.8%7 ≤ 𝑧 ≤ SegU-Net shows better precision in recovering shapes of the ion-ized regions compared to Super-Pixel method. As expected, most ofthe network uncertainty is located on the boundary of neutral regionsor between two large ionized bubbles, when these are percolating andthe gap is getting narrower. This uncertainty has a direct bearing onsmall neutral islands, of a few Mpc scale, residing in vast ionizedregions.
To compare the predicted ionized fields from the 21-cm images math-ematically, we use the Matthews correlation coefficient (MCC) (alsoknown as 𝑟 𝜙 coefficient) defined as: 𝑟 𝜙 = 𝑁 TP · 𝑁 TN − 𝑁 FP · 𝑁 FN √︁ ( 𝑁 TP + 𝑁 FP )( 𝑁 TP + 𝑁 FN )( 𝑁 TN + 𝑁 FP )( 𝑁 TN + 𝑁 FN ) , (4)where 𝑁 TP and 𝑁 FN are the total number of neutral and ionized pixelsrecovered correctly respectively. 𝑁 TN is the total number of pixelsincorrectly guessed as neutral and 𝑁 FP is the total number of pixelsincorrectly guessed as ionized. MCC is a useful metric to correlatebinary fields.In Figure 4 we show the MCC estimated from the fields segmentedinto ionized and neutral regions in our testing sets. In the left panel MNRAS , 1–15 (2021) Figure 5.
Left panel:
Comparison of the simulated neutral fraction against the recovered one. Error-bar and color-bar are the same as in the previous figure.
Right panel:
The same comparison for the ‘fiducial’ simulation. We also include the results from C RAY simulation. The redshift of C RAY simulation are 𝑧 = . , . , . , .
40 corresponding to a volume averaged neutral fraction of 𝑥 vHI = . , . , . , .
70. The green dots with relative confidenceinterval are predictions performed with
SegU-Net and green squares are with the Super-Pixel method. we provide a scatter plot of MCC values against the reionizationhistory ( 𝑥 vHI ) for the ‘random’ testing set. We indicate the redshift ofthe realization by the color of the points and respective confidenceinterval with an error-bar. In an inset panel, we show the number ofsamples in our training set at different neutral fraction. After a firstattempt, we realized that to overcome the class unbalanced problemwe required a larger representation of simulations at early ( 𝑥 vHI ≈ 𝑥 vHI ≈ 𝑥 vHI ≈ . . 𝑟 𝜙 value for the overall prediction data set (Figure 4,left panel) is about 87 per cent for SegU-Net (blue dashed line) and62 per cent in the case of the Super-Pixel method (orange dashedline). The noise level increases with redshift, therefore the score isslightly less accurate for redshift 𝑧 ≥ .
25 with accuracy 85 per cent,meanwhile higher for lower redshift 𝑧 ≤ .
75 with 88 per cent. Inthe future, we consider to increase the proportion of the training datawith high redshift to decrease this performance dissimilarity. Thesame trend is present in the case of the Super-Pixel method, withaccuracy of 60 per cent and 63 per cent respectively.In the right panel of Figure 4, we compare the MCC values from
SegU-Net (blue line with circles) with that from the Super-Pixelmethod (orange line with squares) for our ‘fiducial’ simulation. Aswe already know from Giri et al. (2018b), the Super-Pixel methodperforms best for 𝑥 VHI ≈ . SegU-Net we are able to overcome this problemby employing a specifically designed BCE loss function during thevalidation process after each training epoch. Therefore, the average 𝑟 𝜙 value for the ‘fiducial’ simulation is about 91 per cent for SegU-Net (blue dashed line) and 70 per cent in the case of the Super-Pixelmethod (orange dashed line). In Table 3, we summarise the 𝑟 𝜙 scorefor the two test sets. After identifying the ionized regions we can determine the volume av-eraged neutral fraction 𝑥 vHI which quantifies the reionization history.In Figure 5 we show the volume averaged neutral fraction 𝑥 vHI , predicted as calculated from the recovered binary fields extracted by the twomethods. In the left panel we show the 𝑥 vHI , predicted of the SegU-Net outputs at true volume averaged neutral fraction 𝑥 vHI , true for our‘random’ testing set. The color of the points indicates the redshifts.The black dashed line indicates 𝑥 vHI , predicted = 𝑥 vHI , true . Except for afew points, all the points lie on or near the black dashed line.In the right panel of Figure 5 we compare the results of 𝑥 VHI , predicted derived with the Super-Pixel method (orange line with squares) and SegU-Net (blue line with circles) for our ‘fiducial’ simulation. Again,the black dashed line represent 𝑥 VHI , predicted = 𝑥 VHI , true . In the caseof our neural network all results lie within the half standard devia-tion (0 . 𝜎 ) of the true value (gray dashed lines). With Super-Pixelmethod, this is true only for 𝑥 VHI ≈ .
6. The recovered neutral frac-tion is either underestimated at 𝑥 VHI > . 𝑥 VHI < . From the 3D tomographic data that will be produced with the upcom-ing SKA experiment, we will be able to study the size distribution ofneutral or ionized region during the EoR. Ionized regions are oftencalled bubbles and whereas neutral regions are referred to as islands.The Bubble and Island size distributions (BSDs and ISDs) are usefulto derive the properties of reionization and its evolution (Xu et al.2017; Giri et al. 2019a). Several approaches ware presented to cal-culate this distribution (Lin et al. 2016; Kakiichi et al. 2017; Giriet al. 2018a). In this work, we employ the Mean-Free-Path method(MFP; Mesinger & Furlanetto 2007; Giri et al. 2018a) to calculatethe size distribution ( 𝑅 d 𝑁 d 𝑅 ) of recovered neutral (ISD) and ionizedfield (BSD). MNRAS000
6. The recovered neutral frac-tion is either underestimated at 𝑥 VHI > . 𝑥 VHI < . From the 3D tomographic data that will be produced with the upcom-ing SKA experiment, we will be able to study the size distribution ofneutral or ionized region during the EoR. Ionized regions are oftencalled bubbles and whereas neutral regions are referred to as islands.The Bubble and Island size distributions (BSDs and ISDs) are usefulto derive the properties of reionization and its evolution (Xu et al.2017; Giri et al. 2019a). Several approaches ware presented to cal-culate this distribution (Lin et al. 2016; Kakiichi et al. 2017; Giriet al. 2018a). In this work, we employ the Mean-Free-Path method(MFP; Mesinger & Furlanetto 2007; Giri et al. 2018a) to calculatethe size distribution ( 𝑅 d 𝑁 d 𝑅 ) of recovered neutral (ISD) and ionizedfield (BSD). MNRAS000 , 1–15 (2021)
M. Bianco et al.
Figure 6.
The size distribution of neutral regions (ISD; left column) the size distribution of ionized region (BSD; right column). Rows from top to bottomrepresents early ( 𝑥 vHI = . 𝑥 vHI = .
5) and late ( 𝑥 vHI = .
2) stages of reionization respectively. On each panel, we show the size distributions fromthe binary fields of ‘fiducial’ simulation recovered by
SegU-Net (blue line) and its respective confidence interval (blue shadow). The size distributions of theground truth and binary field recovered by the Super-Pixel method are given by black dashed lines and orange lines. On the bottom of each size distributionpanel, we show the relative deviation from the true binary field distribution.MNRAS , 1–15 (2021) Figure 7.
Dimensionless power spectra of the reference simulation recoveredneutral field by our network (blue line) and its respective confidence interval(blue shadow). Compared at early, middle and late stage of reionization (fromtop to bottom 𝑥 VHI = . , . , .
2) with the same quantity derived from theground truth (black dashed line) and the Super-Pixel method (orange line).On the bottom of each panel, we show the relative difference compared to theground truth for both cases, network and Super-Pixel method.
Figure 8.
Comparison of topology of the identified regions with Betti num-bers estimated from the original neutral field (black dashed line), the
SegU-Net (blue line with circles) and the Super-Pixel method (orange line withsquares), for the case of our ‘fiducial’ simulation. The top, middle and bot-tom panels give 𝛽 , 𝛽 and 𝛽 respectively. The Betti numbers recoveredwith SegU-Net matches with the ground truth better compared to the onesrecovered with the Super-Pixel method. MNRAS000
SegU-Net (blue line with circles) and the Super-Pixel method (orange line withsquares), for the case of our ‘fiducial’ simulation. The top, middle and bot-tom panels give 𝛽 , 𝛽 and 𝛽 respectively. The Betti numbers recoveredwith SegU-Net matches with the ground truth better compared to the onesrecovered with the Super-Pixel method. MNRAS000 , 1–15 (2021) M. Bianco et al.
Figure 9.
On the left, the Matthews correlation coefficient 𝑟 𝜙 of the recovered binary field against its volume averaged neutral fraction. We compare theprediction set for a high noise level ( 𝑡 obs = ℎ , green line with squares) and low noise level ( 𝑡 obs = ℎ , red line with triangles) against the noise levelemployed during the training ( 𝑡 obs = ℎ , blue line with squares). Horizontal dashed lines of the respective color represent the MCC average score of thereference simulation. On the right, the evolution of the MCC score for increasing observation time for a set of mock observations with volume averaged neutralfraction of 𝑥 VHI = . ( 𝑧 = . ) , . ( 𝑧 = . ) and 0 . ( 𝑧 = . ) , respectively in blue, orange and green color. In the same panel, an inset plot showsthe signal-to-noise ratio (SNR = 𝜎 / 𝜎 noise ) of 21-cm images at resolution corresponding to a maximum baseline of 2 km we achieve for different observationtimes. In Figure 6 left column, we show the ISDs and on the right col-umn we show the BSDs of the binary fields recovered with
SegU-Net (blue line) and Super-Pixel method (orange line) compared to theground truth (black dashed line). The Super-Pixel method performsbest when simulation are halfway through the reionization process 𝑥 vHI = . SegU-Net . We show difference percentage relative tothe ground truth in plots below the ISDs. The blue shaded regionshows the error on each of the ISDs determined by
SegU-Net . Inboth the ISD and BSD case, the main difference between the tworecovered distribution occurs at the earlier 𝑥 vHI = . 𝑥 vHI = . SegU-Net shows a rela-tive difference of a few percentage while the distributions determinedfrom the Super-Pixel segmentations show relative differences up to10 per cent for large sizes.
The dimensionless power spectrum of the neutral field is definedas Δ
2x x = 𝑘 𝑃 x x ( 𝑘 )/ 𝜋 , where 𝑃 x x is the auto-power spectrumthat quantifies the fluctuations due to the distribution of neutral re-gions. These fluctuations contribute to the 21-cm power spectrumthat is observed with radio interferometric telescopes. See for exam-ple Furlanetto et al. (2006) and Lidz et al. (2007) for descriptions ofthe fluctuations of the 21-cm signal. In this section, we study the Δ
2x x estimated from the neutral fields recovered from various methods.In Figure 7, we consider the ‘fiducial’ simulation at three stagesof reionization, which are 𝑥 VHI = . 𝑥 VHI = . 𝑥 VHI = . 𝑘 < . − with relative difference within 25 per cent for lowerk-values. The Δ
2x x of the neutral field recovered by the Super-Pixelmethod at early and late times have the correct shape but differs inmagnitude. The Δ
2x x of the neutral field recovered by
SegU-Net per-form quite well at all three stages of reionization. The network main- tains a maximum difference, when compared with the ground truth,of a few tens of per cent at all scales. For 𝑘 (cid:46) . − , the networkuncertainty interval grows to 25-50 per cent relative difference. During reionization ionized bubbles form, grow and connect witheach other to form a complex topology (Furlanetto & Oh 2016).Various studies have proposed topological descriptors for this distri-bution, such as Euler characteristics (e.g. Friedrich et al. 2011) andBetti numbers (Giri & Mellema 2020; Kapahtia et al. 2021). Giri &Mellema (2020) found that Betti numbers contain more informationcompared to the Euler characteristics. Therefore in this section westudy the zeroth 𝛽 , first 𝛽 and second 𝛽 Betti number (Betti 1870)of the binary 3D maps recovered by various feature identificationmethods. 𝛽 , 𝛽 and 𝛽 describe the number of isolated ionizing regions,tunnels and isolated neutral regions respectively. In the top, middleand bottom panels of Figure 8, we show the 𝛽 , 𝛽 and 𝛽 valuesestimated from the recovered binary fields of our ‘fiducial’ model at 𝑥 VHI between 0.1 and 0.9. The black, blue and orange curves representthe Betti numbers calculated from the ground truth, recovered fieldwith
SegU-Net and Super-Pixel method respectively. In line withthe results for the other quantities discussed above, we find that thetopology recovered with
SegU-Net is much closer to the groundtruth than the one recovered by the Super-Pixel method.
We have trained and tested
SegU-Net for one specific noise level,corresponding to the theoretically expected noise for 𝑡 obs = MNRAS , 1–15 (2021) the theoretical noise level is not achieved due to complications withthe calibration. It is therefore important to test to which extent theperformance of our network is sensitive to the noise level in the ac-tual data. To change the noise level we choose different observingtimes, one shorter ( 𝑡 obs =
500 h) and one longer ( 𝑡 obs = √ √︁ / 𝑟 𝜙 coefficient of therecovered binary field against the volume averaged neutral fraction 𝑥 VHI . We compare the prediction on the reference simulation for thehigher ( 𝑡 obs =
500 h, green line with squares) and lower noise case( 𝑡 obs = 𝑡 obs = 𝑥 VHI > .
8, due to the redshift dependency of the simulated noise.We also want to test how far we can push our
SegU-Net trainedon data with 𝑡 obs = 𝑟 𝜙 coefficient at different observation times 𝑡 obs ,for three different stages of reionization in our reference simulation,namely for volume averaged neutral fractions 𝑥 VHI is 0.2 (blue linewith squares), 0.5 (orange line with circles) and 0.8 (green line withtriangles), corresponding to redshift 𝑧 = . , .
032 and 8 . 𝑡 obs (cid:38)
500 h,where 𝑟 𝜙 (cid:38) .
85. The spike in the curve for 𝑥 VHI = . 𝑡 obs = = 𝜎 / 𝜎 noise (e.g. Kakiichi et al. 2017), where = 𝜎 and 𝜎 noise are the standarddeviations of the 21-cm signal and noise respectively at the resolu-tion corresponding to a maximum baseline of 2 km. From this weconclude that a good performance ( 𝑟 𝜙 (cid:38) .
85) requires a SNR (cid:38) To train
SegU-Net , we relied on for creating the train-ing set. However, our Universe may not exactly be described by thissemi-numerical code. If our neural network has learnt to find struc-tures in simulations only, then we cannot use it for SKAobservations. To make sure that the neural network is not over-fitted,we consider a different reionization simulation code to built the mockobservations. We describe this simulation in § 6.1 and present ourresults in § 6.2.
We first simulate the matter density field and track the evolutionof cosmic structures by using the
CUBEP M 𝑁 -body code (Harnois-Déraps et al. 2013). The simulation is carried out in a volume of ( ) with 64 billion particles. Dark matter haloes down toa mass of 10 𝑀 (cid:12) is found at various redshift using the sphericalaverage halo-finder (Watson et al. 2013), meanwhile haloes with Figure 10.
An example of the smoothed cube slice from C RAY simulationon the right, employed to test the stability and reliability of the network. At 𝑧 = .
06 corresponding to a volume averaged neutral fraction of 𝑥 VHI = . mass between 10 and 10 𝑀 (cid:12) are implemented with a sub-gridmethod (Ahn et al. 2015). We use the same cosmology that is givenin § 2.1.We then employ the C RAY radiative transfer (RT) code (Mellemaet al. 2006b) to simulate the cosmic reionization. C RAY requires thematter density field in 3D grid. Therefore the distribution 𝑁 -bodyparticles are put in 3D grids with a smoothed particle hydrodynamicmethod (e.g. Shapiro et al. 1996; Mao et al. 2019). This grid hasspatial resolution of Δ 𝑥 = . mesh-grid. Allhaloes emits 4 × s − photons per unit time. We consider partialsuppression only on haloes with mass between 10 and 10 𝑀 (cid:12) .When these are located in a ionized cell (above 10%) they emits10 s − photons per unit time. C RAY outputs the hydrogen ionizationfield at a time interval of 11.5 Myrs. For more details on the RT and 𝑁 -body simulations methods, see Iliev et al. (2012) and Bianco et al.(2021).We derived the differential brightness temperature 𝛿𝑇 b from theionization field and the density using Equation 2. We selectedfour outputs, which are at redshifts 𝑧 = . , . , . , . 𝑥 VHI = . , . , . , .
70 respectively. Simulated 𝛿𝑇 b from these epochswere converted into mock observations (see § 2.3). We use thesemock observations as a testing set.In Figure 10 right panel, we show a slice of the calculated 𝛿𝑇 b forredshift 𝑧 = .
06 ( 𝑥 VHI = . 𝐵 = 𝛿𝑇 b data set. We applied our network to mock 𝛿𝑇 𝑏 cubes calculated with the C RAY code, with a spatial resolution close to the simula-tions employed in the training process, ∼ MNRAS000
06 ( 𝑥 VHI = . 𝐵 = 𝛿𝑇 b data set. We applied our network to mock 𝛿𝑇 𝑏 cubes calculated with the C RAY code, with a spatial resolution close to the simula-tions employed in the training process, ∼ MNRAS000 , 1–15 (2021) M. Bianco et al. ered binary field, similar to the results in § 4.1, is shown in Figure 10.In the left panel, the red/blue color indicates the network predictionWe go through the same process presented in § 4. The 𝑟 𝜙 score isrepresented by the green error-bar points on the right panel of Fig-ure 4, from left to right we have redshift 𝑧 = . , . , . , . 𝑥 VHI = . , . , . , .
70, green squares represent the scoreobtained with the Super-Pixel method. As we can see our neuralnetwork is performing with the same accuracy as with the semi-numerical simulation employed during the validation and training,discussed in § 4.1. For 𝑥 VHI ≈ . SegU-Net performs slightly betterthan the Super-Pixel method. Super-Pixel methods shows a drop inaccuracy at the late stages of reionization ( 𝑥 VHI < . 𝑥 VHI ,in Figure 5 right panel, the green error-bar points are the same dataas mentioned above. As we can see, also for the C RAY simulation therecovered quantity resides within the 0.5- 𝜎 confidence interval. Forthe Super-Pixel results this is true only for its best performing case,when 𝑥 VHI ≈ .
5, where it achieves approximately the same precisionas
SegU-Net . In this work we have developed a convolutional neural network basedon the U-Net architecture which can be used to segment redshifted21-cm image observations into neutral and ionized regions. We haveshown that this application of deep learning provides a fast andstable method that significantly improves the identification of ion-ized/neutral regions during the epoch of reionization over previouslyproposed methods. To train our network, we employ a large set ofsimulated mock observations of the 21-cm signal.Our image segmentation network,
SegU-Net , also contains anuncertainty estimator. This uncertainty estimator is a simple but effi-cient application of the test-time augmented (TTA) technique. Withthis uncertainty estimator, our network can create a pixel by pixelerror map during the prediction process. The pixel by pixel errormap can later be used to determine the error in any quantity derivedfrom the segmentation.We compare the accuracy of our approach with a feature findingmethod from the field of image processing, the so-called Super-Pixelmethod which (Giri et al. 2018b) showed to be superior over simplethresholding methods. The results show that our neural network canidentify neutral regions in the mock observations at least as well andoften much better than the Super-Pixel method. We show a visualcomparison and the resulting pixel per pixel error map tested on our‘fiducial’ model. This error map gives an useful insight about theparts of the image that are hard to recover and also helps us checkfor over-fitting.We studied the accuracy of a range of derived quantities fromthe recovered binary fields, comparing the performance of
SegU-Net with the Super-Pixel method. These quantities are the volumeaveraged ionization fraction — the evolution of which provides thereionization history, the size distribution of the ionized (BSD) andneutral (ISD) regions, the dimensionless power spectra of the re-covered binary fields and the three Betti numbers, which quantifythe topology of the segmented data sets. For all quantities we findthat the
SegU-Net results are more accurate than the Super-Pixelresults, especially for the early and late stages of reionization wherethe Super-Pixel method often struggles to produce accurate results.Machine learning methods generally are sensitive to the propertiesof the training set. We therefore tested
SegU-Net on input data with different properties than the training set. First, we analysed theperformance on data sets with different noise levels than the networkwas trained for and found that
SegU-Net yields accurate results fordata sets in which the noise level is characterised by an observing timeof 𝑡 obs > ℎ , which approximately corresponds to a SNR (cid:38) SegU-Net achieves the same level of accuracy when appliedto this data set and therefore is not sensitive to the type of simulationemployed during the training process.We want to point out that similar efforts are being made by Gagnon-Hartman et al (in prep). They focus on reconstructing the segmentedmaps of ionized and neutral regions in the context of foreground mit-igation using the foreground avoidance method (e.g. Liu & Tegmark2011; Pober et al. 2014), and also consider the possibility of do-ing so with instruments that are not optimized for imaging such asHERA. We include the effect of instrumental noise and study in greatdetail the summary statistics of the reconstructed binary maps andthe dependency of the results on the noise level. In the future, wewill extend our study to include the impact of foreground mitigationstrategies while recovering the summary statistics.Here we assumed the spin temperature to be saturated ( 𝑇 S (cid:29) 𝑇 CMB ). However, it is possible to have such a scenario where thisassumption fail, specially during the time when reionization starts.In the furture, we will evolve our
SegU-Net to identify H ii regionsin such scenarios. Even though our network is built to identify H iiregions, U-Net architecture can be trained to identify any interestingfeature. Before reionization started, the luminous sources heated theIGM and leave its impact on the 21-cm signal (e.g. Ross et al. 2017,2019). The U-Net architecture can also be trained to identify theseheated regions.
Tensorflow (Abadi et al. 2015) and
Keras (Chollet et al. 2017). The algorithms and image processing toolsoperated on our data were performed with the help of
NumPy (Har-ris et al. 2020),
SciPy (Virtanen et al. 2020) and scikit-image packages (van der Walt et al. 2014). All figures were created with mathplotlib (Hunter 2007).
MNRAS , 1–15 (2021) Data Availability:
The data underlying this article is available uponrequest, and can also be re-generated from scratch using the pub-licly available , CUBEP M , C RAY and
Tools21cm code.The
SegU-Net code and its trained network weights are avail-able on the author’s
GitHub page: https://github.com/micbia/SegU-Net . REFERENCES
Abadi M., et al., 2015, TensorFlow, http://tensorflow.org/
Ahn K., Iliev I. T., Shapiro P. R., Srisawat C., 2015, Monthly Notices of theRoyal Astronomical Society, 450, 1486Betti E., 1870, Annali di Matematica Pura ed Applicata (1867-1897), 4, 140Bianco M., Iliev I. T., Ahn K., Giri S. K., Mao Y., Park H., Shapiro P. R.,2021 ( arXiv:2101.01712 )Chardin J., Uhlrich G., Aubert D., Deparis N., Gillet N., Ocvirk P., Lewis J.,2019, MNRAS, 490, 1055–1065Charnock T., Perreault-Levasseur L., Lanusse F., 2020, Bayesian Neural Net-works ( arXiv:2006.01490 )Chen B. H., et al., 2020, MNRAS, 501, 3951–3961Chollet F., Allaire J., et al., 2017, Keras, https://github.com/rstudio/keras
Cohen A., Fialkov A., Barkana R., Monsalve R., 2019, arXiv e-prints, p.arXiv:1910.06274Cui Y., Jia M., Lin T., Song Y., Belongie S. J., 2019 ( arXiv:1901.05555 )Datta K. K., Bharadwaj S., Choudhury T. R., 2007, Monthly Notices of theRoyal Astronomical Society, 382, 809Davies F. B., et al., 2018, The Astrophysical Journal, 864, 142DeBoer D. R., et al., 2017, Publications of the Astronomical Society of thePacific, 129, 45001Dumoulin V., Visin F., 2016 ( arXiv:1603.07285 )Fan X., et al., 2006, The Astronomical Journal, 132, 117–136Ferrara A., Pandolfi S., 2014, Reionization of the Intergalactic Medium( arXiv:1409.4946 )Friedrich M. M., Mellema G., Alvarez M. A., Shapiro P. R., Iliev I. T., 2011,MNRAS, 413, 1353Furlanetto S. R., Oh S. P., 2016, Monthly Notices of the Royal AstronomicalSociety, 457, 1813Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1–15Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Physics Reports, 433, 181Ghara R., Choudhury T. R., Datta K. K., Choudhuri S., 2016, MNRAS, 464,2234–2248Ghara R., et al., 2020, Monthly Notices of the Royal Astronomical Society,493, 4728Gillet N., Mesinger A., Greig B., Liu A., Ucci G., 2019, MNRASGiri S. K., 2019, PhD thesis, Department of Astronomy, Stockholm UniversityGiri S. K., Mellema G., 2020, arXiv e-prints, p. arXiv:2012.12908Giri S. K., Mellema G., Dixon K. L., Iliev I. T., 2018a, Monthly Notices ofthe Royal Astronomical Society, 473, 2949Giri S. K., Mellema G., Ghara R., 2018b, MNRAS, 479, 5596–5611Giri S. K., Mellema G., Aldheimer T., Dixon K. L., Iliev I. T., 2019a, MonthlyNotices of the Royal Astronomical Society, 489, 1590Giri S. K., D’Aloisio A., Mellema G., Komatsu E., Ghara R., Majumdar S.,2019b, Journal of Cosmology and Astroparticle Physics, 2019, 058Giri S. K., Mellema G., Jensen H., 2020a, Journal of Open Source Software,5, 2363Giri S. K., Zackrisson E., Binggeli C., Pelckmans K., Cubo R., 2020b, MN-RAS, 491, 5277Glorot X., Bordes A., Bengio Y., 2011, in Gordon G., Dunson D., Dudík M.,eds, Proceedings of Machine Learning Research Vol. 15, Proceedingsof the Fourteenth International Conference on Artificial Intelligence andStatistics. JMLR Workshop and Conference Proceedings, Fort Laud-erdale, FL, USA, pp 315–323, http://proceedings.mlr.press/v15/glorot11a.html
Greig B., Mesinger A., Bañados E., 2019, MNRAS, 484, 5094Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., StuartJ., Wyithe B., 2020a, arXiv:2008.02639 Greig B., et al., 2020b, arXiv:2006.03203, 000, 0Guzman E., Meyers J., 2021 ( arXiv:2101.01214 )Harnois-Déraps J., Pen U. L., Iliev I. T., Merz H., Emberson J. D., DesjacquesV., 2013, MNRAS, 436, 540Harris C. R., et al., 2020, Nature, 585, 357–362Hassan S., Andrianomena S., Doughty C., 2020, MNRAS, 494, 5761–5774Hinton G. E., Srivastava N., Krizhevsky A., Sutskever I., Salakhutdinov R. R.,2012 ( arXiv:1207.0580 )Hortúa H. J., Volpi R., Malagò L., 2020 ( arXiv:2005.02299 )Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Ichikawa K., Barkana R., Iliev I. T., Mellema G., Shapiro P. R., 2010, MN-RAS, 406, 2521Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., AhnK., 2012, Monthly Notices of the Royal Astronomical Society, 423,2222–2253Ioffe S., Szegedy C., 2015 ( arXiv:1502.03167 )Jarrett K., Kavukcuoglu K., Ranzato M., LeCun Y., 2009, in 2009 IEEE12th International Conference on Computer Vision. pp 2146–2153,doi:10.1109/ICCV.2009.5459469Jeffrey N., Lanusse F., Lahav O., Starck J.-L., 2020, MNRAS, 492, 5023–5029Jennings W. D., Watkinson C. A., Abdalla F. B., McEwen J. D., 2018, MN-RAS, 483, 2907–2922Kakiichi K., et al., 2017, Monthly Notices of the Royal Astronomical Society,471, 1936Kapahtia A., Chingangbam P., Ghara R., Appleby S., Choudhury T. R., 2021( arXiv:2101.03962 )Kern N. S., Liu A., Parsons A. R., Mesinger A., Greig B., 2017, The Astro-physical Journal, 848, 23Kingma D. P., Ba J., 2014 ( arXiv:1412.6980 )Komatsu E., et al., 2009, The Astrophysical Journal Supplement Series, 180,330–376Komatsu E., et al., 2011, The Astrophysical Journal Supplement Series, 192,18Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilo-metre Array (AASKA14). http://pos.sissa.it/
Lee C., 2019, Planetary and Space Science, 170, 16–28Li W., et al., 2019, MNRAS, 485, 2628Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007,The Astrophysical Journal, 659, 865Lin Y., Oh S. P., Furlanetto S. R., Sutter P. M., 2016, Monthly Notices of theRoyal Astronomical Society, 461, 3361Liu A., Tegmark M., 2011, Phys. Rev. D, 83, 103006Long J., Shelhamer E., Darrell T., 2014 ( arXiv:1411.4038 )Madau P., Meiksin A., Rees M. J., 1997, The Astrophysical Journal, 475, 429Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A., Bharadwaj S.,Mellema G., 2018, Monthly Notices of the Royal Astronomical Society,476, 4007Makinen T. L., Lancaster L., Villaescusa-Navarro F., Melchior P., Ho S.,Perreault-Levasseur L., Spergel D. N., 2020 ( arXiv:2010.15843 )Mangena T., Hassan S., Santos M. G., 2020, MNRAS, 494, 600–606Mao Y., Koda J., Shapiro P. R., Iliev I. T., Mellema G., Park H., Ahn K.,Bianco M., 2019, MNRAS, 22, 1Mason C. A., Gronke M., 2020, MNRAS, 499, 1395McGreer I. D., Mesinger A., Fan X., 2011, MNRAS, 415, 3237–3246McGreer I. D., Mesinger A., D’Odorico V., 2014, MNRAS, 447, 499–505Mehta P., Bukov M., Wang C.-H., Day A. G., Richardson C., Fisher C. K.,Schwab D. J., 2019, Physics Reports, 810, 1–124Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006a, MNRAS, 372, 679Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006b, Monthly Notices ofthe Royal Astronomical Society, 372, 679Mellema G., et al., 2013, Experimental Astronomy, 36, 235Mertens F. G., et al., 2020, Monthly Notices of the Royal AstronomicalSociety, 493, 1662Mesinger A., Furlanetto S., 2007, The Astrophysical Journal, 669, 663–675Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955Mondal R., et al., 2020, Monthly Notices of the Royal Astronomical Society,498, 178 MNRAS000
Lee C., 2019, Planetary and Space Science, 170, 16–28Li W., et al., 2019, MNRAS, 485, 2628Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007,The Astrophysical Journal, 659, 865Lin Y., Oh S. P., Furlanetto S. R., Sutter P. M., 2016, Monthly Notices of theRoyal Astronomical Society, 461, 3361Liu A., Tegmark M., 2011, Phys. Rev. D, 83, 103006Long J., Shelhamer E., Darrell T., 2014 ( arXiv:1411.4038 )Madau P., Meiksin A., Rees M. J., 1997, The Astrophysical Journal, 475, 429Majumdar S., Pritchard J. R., Mondal R., Watkinson C. A., Bharadwaj S.,Mellema G., 2018, Monthly Notices of the Royal Astronomical Society,476, 4007Makinen T. L., Lancaster L., Villaescusa-Navarro F., Melchior P., Ho S.,Perreault-Levasseur L., Spergel D. N., 2020 ( arXiv:2010.15843 )Mangena T., Hassan S., Santos M. G., 2020, MNRAS, 494, 600–606Mao Y., Koda J., Shapiro P. R., Iliev I. T., Mellema G., Park H., Ahn K.,Bianco M., 2019, MNRAS, 22, 1Mason C. A., Gronke M., 2020, MNRAS, 499, 1395McGreer I. D., Mesinger A., Fan X., 2011, MNRAS, 415, 3237–3246McGreer I. D., Mesinger A., D’Odorico V., 2014, MNRAS, 447, 499–505Mehta P., Bukov M., Wang C.-H., Day A. G., Richardson C., Fisher C. K.,Schwab D. J., 2019, Physics Reports, 810, 1–124Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006a, MNRAS, 372, 679Mellema G., Iliev I. T., Pen U.-L., Shapiro P. R., 2006b, Monthly Notices ofthe Royal Astronomical Society, 372, 679Mellema G., et al., 2013, Experimental Astronomy, 36, 235Mertens F. G., et al., 2020, Monthly Notices of the Royal AstronomicalSociety, 493, 1662Mesinger A., Furlanetto S., 2007, The Astrophysical Journal, 669, 663–675Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955Mondal R., et al., 2020, Monthly Notices of the Royal Astronomical Society,498, 178 MNRAS000 , 1–15 (2021) M. Bianco et al.
Murray S. G., Greig B., Mesinger A., Muñoz J. B., Qin Y., Park J., WatkinsonC. A., 2020, Journal of Open Source Software, 5, 2582Ota K., et al., 2008, The Astrophysical Journal, 722, 803Ouchi M., et al., 2010, The Astrophysical Journal, 723, 869Perez L., Wang J., 2017 ( arXiv:1712.04621 )Planck Collaboration 2019, A&APlanck Collaboration et al., 2016, preprintPober J. C., et al., 2014, \apj, 782, 66Pritchard J. R., Furlanetto S. R., 2007, Monthly Notices of the Royal Astro-nomical Society, 376, 1680Robertson B. E., Ellis R. S., Furlanetto S. R., Dunlop J. S., 2015, The Astro-physical Journal, 802, L19Ronneberger O., Fischer P., Brox T., 2015 ( arXiv:1505.04597 )Ross H. E., Dixon K. L., Iliev I. T., Mellema G., 2017, Monthly Notices ofthe Royal Astronomical Society, 468, 3785Ross H. E., Dixon K. L., Ghara R., Iliev I. T., Mellema G., 2019, MNRAS,487, 1101Sadr A. V., Farsian F., 2020 ( arXiv:2004.04177 )Salehi S. S. M., Erdogmus D., Gholipour A., 2017 ( arXiv:1706.05721 )Schmit C. J., Pritchard J. R., 2018, \mnras, 475, 1213Schroeder J., Mesinger A., Haiman Z., 2013, MNRAS, 428, 3058Shapiro P. R., Martel H., Villumsen J. V., Owen J. M., 1996, ApJ, pp 270–330Simonyan K., Zisserman A., 2015 ( arXiv:1409.1556 )Srivastava N., Hinton G., Krizhevsky A., Sutskever I., Salakhutdinov R.,2014, Journal of Machine Learning Research, 15, 1929Sullivan D., Iliev I. T., Dixon K. L., 2017, MNRAS, 473, 38–58Szegedy C., Vanhoucke V., Ioffe S., Shlens J., Wojna Z., 2015( arXiv:1512.00567 )Tingay S. J., et al., 2013, Publications of the Astronomical Society of Australia(PASA), 30, 7Totani T., Aoki K., Hattori T., Kawai N., 2016, Publications of the Astronom-ical Society of Japan, 68, 1Trott C. M., et al., 2020, Monthly Notices of the Royal Astronomical Society,493, 4711Villanueva-Domingo P., Villaescusa-Navarro F., 2020 ( arXiv:2006.14305 )Virtanen P., et al., 2020, Nature Methods, 17, 261Wang Y., Huang G., Song S., Pan X., Xia Y., Wu C., 2020( arXiv:2007.10538 )Watkinson C. A., Pritchard J. R., 2015, MNRAS, 454, 1416Watson W. A., Iliev I. T., D’Aloisio A., Knebe A., Shapiro P. R., Yepes G.,2013, Monthly Notices of the Royal Astronomical Society, 433, 1230Wyithe S., Geil P. M., Kim H., 2015, Proceedings of Advancing As-trophysics with the Square Kilometre Array (AASKA14). 9 -13June, 2014. Giardini Naxos, Italy. Online at http://pos.sissa.it/cgi-bin/reader/conf.cgi?confid=215, id.15Xu Y., Yue B., Chen X., 2017, The Astrophysical Journal, 844, 117Yoshiura S., Shimabukuro H., Hasegawa K., Takahashi K., 2020( arXiv:2004.09206 )Zackrisson E., et al., 2020, Monthly Notices of the Royal Astronomical So-ciety, 493, 855Zaroubi S., 2012, Astrophysics and Space Science Library, p. 45–101Zeiler M. D., Fergus R., 2013 ( arXiv:1311.2901 )Zel’Dovich Y. B., 1970, A&A, 500, 13van Haarlem M. P., et al., 2013, Astronomy & Astrophysics, 556, A2van der Walt S., et al., 2014, PeerJ, 2, e453
APPENDIX A:
SegU-Net
HIDDEN LAYER OUTPUTS
We test
SegU-Net to see if it can recover the binary field for asimple case, namely a single spherical neutral region. We assumea uniform density field at 𝑧 = .
032 and calculate the differentialbrightness temperature with Equation 2, adding noise correspondingto 𝑡 obs = 𝐵 = SegU-Net (right panel). The black contour in the left
Figure A1.
Test
SegU-Net on a spherical ionized region.
Left panel : slicethrough the input image. The color map shows the differential brightnesstemperature and the black contour shows the boundary between neutral andionized regions.
Right panel : the recovered binary field with
SegU-Net . Thegreen contour represents again the same boundary. The red region indicatesthe identified neutral regions and blue as ionized region.
Figure A2.
Visual representation of
SegU-Net ’s low dimensional latentspace (bottom layer), which contains information about the extracted featuresof our test input image. panel and green contour in the right panel show the true boundaryof the sphere. For this test
SegU-Net achieves an accuracy of 98percent.In Figure A2 we show the output of the bottom hidden layer of
SegU-Net , which is the last layer of the left part of the U-shapedin the Figure 2. The colour coding is such that blue correspond tonegative, red to positive and white to zero values. This output givesa visual representation of the low dimensional latent space of ournetwork encoder. In our case, this consists of 256 images, whereeach corresponds to a convolutional kernel and contains informationabout the image’s extracted features. The encoder path contracts theinput image from an original mesh size of 128 down to 8 . Theinformation contained in the latent space is then expanded by SegU-Net decoder into a binary map of the same size as the input image(see right panel of Figure A1).
MNRAS , 1–15 (2021) Figure A3.
Example of skip connection between encoder and decoder levels. The top panel shows the architecture of our network. The bottom panels displaythe output of three hidden layers. On the left (green dashed line), a convolutional block (
ConvBlock ) output is interconnected with the output of the second tolast up-sampling operation (central panel, red dashed line). The right most panel shows the result of the merge after a convolution block (black dashed line).
APPENDIX B: SKIP CONNECTION BETWEEN ENCODERAND DECODER LEVELS
The main advantage of a U-shaped network is that it overcomesthe bottleneck limitation encountered by auto-encoder networks (aclassical encoder/decoder architecture), by adding interconnectionsbetween the descending (encoder) and ascending (decoder) paths(Long et al. 2014; Ronneberger et al. 2015). These connections allowfeature representations to pass through the bottleneck (bottom layer)and avoid loss of information due to contraction.In Figure A3, we show a visual example of interconnections be-tween the encoder (left part of the U-shape) and the decoder (rightpart). The top panel shows a schematic representation of our networkarchitecture and the bottom part displays a visual output of threehidden layers, for the test case of a sphere. The left-most panel (connected by a green dashed line) shows theoutput of the second convolutional block in the encoder’s secondlevel. This consists of 32 kernels with a mesh size of 64 . At thislevel the shape and form of the input image are still visible. The centrepanel (connected by a red dashed line) shows the result of the secondto last up-sampling operation of the decoder. The number of kernelsand mesh size match with the corresponding encoder layers. The skipconnection merges the encoder and decoder-level output, for a totalof 64 images with mesh 64 . The right-most panel (connected witha black dashed line) shows the concatenation after a convolutionalblock. The effect of the up-sampling operation is still visible. MNRAS000