Gu J., Xu H., Wang J., An T., Chen W., 2013, ApJ, 773, 38

Harker G. et al., 2009, MNRAS, 397, 1138

Harnois-Déraps J., Pen U.-L., Iliev I. T., Merz H., Emberson J. D., Desjacques V., 2013, MNRAS, 436, 540

Harris C. R. et al., 2020, Nature, 585, 357

Hirling P., Bianco M., Giri S. K., Iliev I. T., Mellema G., Kneib J.-P., 2023, preprint (arXiv:2311.01492)

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Hütsi G., Raidal M., Urrutia J., Vaskonen V., Veermäe H., 2023, Phys. Rev. D, 107, 043502

Hutter A., 2018, MNRAS, 477, 1549

Iliev I. T., Mellema G., Pen U. L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625

Iliev I. T., Mellema G., Shapiro P. R., Pen U.-L., Mao Y., Koda J., Ahn K., 2012, MNRAS, 423, 2222

Jelić V. et al., 2008, MNRAS, 389, 1319

Kapahtia A., Chingangbam P., Appleby S., 2019, J. Cosmol. Astropart. Phys., 2019, 053

Kapahtia A., Chingangbam P., Ghara R., Appleby S., Choudhury T. R., 2021, J. Cosmol. Astropart. Phys., 2021, 026

Kerrigan J. R. et al., 2018, ApJ, 864, 131

Kingma D. P., Ba J., 2014, preprint (arXiv:1412.6980)

Komatsu E. et al., 2009, ApJ, 180, 330

Koopmans L. V. E. et al., 2015, Proc. Sci. The Cosmic Dawn and Epoch of Reionization with the Square Kilometre Array. SISSA, Trieste, PoS(AASKA14)001

Li X., Chen S., Hu X., Yang J., 2018, preprint (arXiv:1801.05134)

Liu A., Shaw J. R., 2020, Publ. Astron. Soc. Pac., 132, 062001

Liu A., Tegmark M., Zaldarriaga M., 2009a, MNRAS, 394, 1575

Liu A., Tegmark M., Bowman J., Hewitt J., Zaldarriaga M., 2009b, MNRAS, 398, 401

Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023019

Madau P., Meiksin A., Rees M. J., 1997, ApJ, 475, 429

Mehra J., Neeru N., 2016, Imperial J. Int. Res., 03, 8

Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2006, New Astron., 11, 374

Mellema G. et al., 2013, Exp. Astron., 36, 235

Mellema G., Koopmans L., Shukla H., Datta K. K., Mesinger A., Majumdar S., 2015, Proc. Sci., HI tomographic imaging of the Cosmic Dawn and Epoch of Reionization with SKA. SISSA, Trieste, PoS(AASKA14)010

Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640

Mertens F. G. et al., 2020, MNRAS, 493, 1662

Mesinger A., Furlanetto S., 2007, ApJ, 669, 663

Mesinger A., Furlanetto S., Cen R., 2011, MNRAS, 411, 955

Morales M. F., Bowman J. D., Cappallo R., Hewitt J. N., Lonsdale C. J., 2006a, New Astron. Rev., 50, 173

Morales M. F., Bowman J. D., Hewitt J. N., 2006b, ApJ, 648, 767

Murray S. G., Trott C. M., 2018, ApJ, 869, 25

Murray S. G., Greig B., Mesinger A., Muñoz J. B., Qin Y., Park J., Watkinson C. A., 2020, J. Open Source Soft., 5, 2582

Naidu R. P. et al., 2022, ApJ, 940, L14

Offringa A. R. et al., 2014, MNRAS, 444, 606

Pawlik A. H., Milosavljević M., Bromm V., 2011, ApJ, 731, 54

Pedregosa F. et al., 2011, J. Mach. Learn. Res., 12, 2825

Perez L., Wang J., 2017, preprint (arXiv:1712.04621)

Planck Collaboration VI, 2020, A&A, 641, A6

Platania P., Bensadoun M., Bersanelli M., De Amici G., Kogut A., Levin S., Maino D., Smoot G. F., 1998, ApJ, 505, 473

Pober J. C. et al., 2014, ApJ, 782, 66

Prelogović D., Mesinger A., Murray S., Fiameni G., Gillet N., 2021, MNRAS, 509, 3852

Pritchard J. R., Furlanetto S. R., 2007, MNRAS, 376, 1680

Ronneberger O., Fischer P., Brox T., 2015, preprint (arXiv:1505.04597)

Ross H. E., Dixon K. L., Iliev' I. T., Mellema G., 2017, MNRAS, 468, 3785

Ross H. E., Dixon K. L., Ghara R., Iliev I. T., Mellema G., 2019, MNRAS, 487, 1101

Ross H. E., Giri S. K., Mellema G., Dixon K. L., Ghara R., Iliev I. T., 2021, MNRAS, 506, 3717

Rumelhart D. E., Zipser D., 1985, Cognitive Sci., 9, 75

Salehi S. S. M., Erdogmus D., Gholipour A., 2017, preprint (arXiv:1706.05721)

Santos M. G., Cooray A., Knox L., 2005, ApJ, 625, 575

Schaeffer T., Giri S. K., Schneider A., 2023, MNRAS, 526, 2942

Schneider A., Schaeffer T., Giri S. K., 2023, Phys. Rev. D, 108, 043030

Selvaraju R. R., Cogswell M., Das A., Vedantam R., Parikh D., Batra D., 2019, Int. J. Comp. Vision, 128, 336

Smirnov O. M., 2011, A&A, 527, A106

Sortino R. et al., 2023, preprint (arXiv:2307.02392)

The HERA Collaboration, 2022a, ApJ, 924, 51

The HERA Collaboration, 2022b, ApJ, 925, 221

Thyagarajan N. et al., 2015, ApJ, 804, 14

Trott C. M. et al., 2020, MNRAS, 493, 4711

Virtanen P. et al., 2020, Nature Methods, 17, 261

van der Walt S. et al., 2014, PeerJ, 2, e453

Wang X., Tegmark M., Santos M. G., Knox L., 2006, ApJ, 650, 529

Wang J. et al., 2013, ApJ, 763, 90

Wang Y., Huang G., Song S., Pan X., Xia Y., Wu C., 2020, preprint (arXiv:2007.10538)

Wang R., Chen Z., Luo Q., Wang F., 2023, preprint (arXiv:2305.09121)

Wyithe S., Geil P. M., Kim H., 2015, Proc. Sci., Imaging HII Regions from Galaxies and Quasars During Reionisation with SKA. SISSA, Trieste, PoS(AASKA14)015

Zackrisson E. et al., 2020, MNRAS, 493, 855

Zaroubi S., 2012, in Wiklind T., Mobasher B., Bromm V., eds, Astrophysics and Space Science Library, Vol 396, The Epoch of Reionization. Springer, Berlin, Heidelberg

## APPENDIX A: HYPERPARAMETER EXPLORATION

As we mentioned in Section 4, we perform an optimization analysis of the SegU-Net hyperparameters. We are aware of tools that automatize the exploratory analysis of the network hyperparameter space, such as Optuna (Akiba et al. 2019). However, constrained by time and computational resources, we manually searched the best-performing parameters through a trial-and-error approach. First, we selected a few combinations of the network parameters and performed a short training of no more than five to ten epochs. Based on the result obtained in this short training, we selected six of the best-performing results with the lowest validation loss and performed a full training to identify the ideal hyperparameters set-up. In future work, we intend to undertake a more comprehensive study.

We list the six best-performing set-ups we tested in Table A1. We include an analysis of seven model parameters, from left to right: the activation function of the convolutional layers, the number of channels for the bottom layer, the number of pooling operations, the dropout rate, the final activation function before the binary operation, the size of the kernel filters and the type of pooling operation. As a loss function, we employed the balanced cross-entropy (BCE) function (Salehi, Erdogmus & Gholipour 2017) and the Adaptive Moment Estimator (Adam) (Kingma & Ba 2014) as the stochastic gradient descent algorithm to minimize the loss. We employed the in bold text for the results presented in this paper. Although we monitored the validation loss to select the best set-up, we noticed that the first and second models gave the worst prediction for images at the edges of the redshift range $z \sim 7$ and 11. The fifth model provided the most balanced result, with an overall $r_\phi \approx 0.7$ score, as shown in Fig. 5 central panel. Moreover, in contrast to the findings by Li et al. (2018), we observe that setting the dropout rate to zero enhances accuracy only for the third-ranked set-up. Meanwhile, other configurations exhibit improved performance when both batch normalization and dropout are included.
<table><caption>**Table A1.** SegU-Net hyperparameter optimization analysis for the best-performing set-ups of seven parameters with optimization on the validation loss.</caption><thead><tr><th>Ranking</th><th>Activation</th><th>Channels latent space</th><th>Depth</th><th>Dropout</th><th>Final activation</th><th>Kernel</th><th>Pooling type</th><th>r<sub>φ</sub> (per cent)</th><th>Validation loss (× 10<sup>-2</sup>)</th></tr></thead><tbody><tr><td>1</td><td>LeakyReLu</td><td>256</td><td>3</td><td>0.42</td><td>σ(x)</td><td>6</td><td>max</td><td>89.08</td><td>6.59</td></tr><tr><td>2</td><td>LeakyReLu</td><td>128</td><td>4</td><td>0.00</td><td>σ(x)</td><td>5</td><td>max</td><td>89.02</td><td>6.62</td></tr><tr><td>3</td><td>Elu</td><td>128</td><td>3</td><td>0.34</td><td>x</td><td>11</td><td>average</td><td>87.76</td><td>7.27</td></tr><tr><td>4</td><td>LeakyReLu</td><td>128</td><td>4</td><td>0.50</td><td>σ(x)</td><td>5</td><td>max</td><td>88.72</td><td>6.85</td></tr><tr><td><b>5</b></td><td><b>ReLu</b></td><td><b>256</b></td><td><b>4</b></td><td><b>0.05</b></td><td><b>x</b></td><td><b>7</b></td><td><b>average</b></td><td><b>88.50</b></td><td><b>7.48</b></td></tr><tr><td>6</td><td>ReLu</td><td>256</td><td>5</td><td>0.14</td><td>x</td><td>7</td><td>max</td><td>86.53</td><td>9.15</td></tr></tbody></table>

## APPENDIX B: INSIDE THE BLACK BOX

The trained model we presented in this paper is able to recover the ionized field from noisy images with residual foreground contamination. This is an indication that the network learns to identify the regions of interest from important hidden features that maximize the recovery. However, the machine learning model’s complexity, high dimensionality and non-linearity make them difficult to interpret and regulate, so these applications are often referred to as a *black box*. Here, we present a first attempt to open and look inside SegU-Net black box. A standard tool to visualize and understand the decisions made by a general convolutional neural network is the Gradient-weighted Class Activation Mapping (Grad-CAM) technique (Selvaraju et al. 2019). This method applied in segmentation and object classification highlights the features of an input image employed in predicting a particular class. Grad-CAM achieves this by computing the gradients of the predicted class $y_c = x_{\text{HI}}^B$ score with respect to the feature maps $F_k$ of all the $k > 0$ convolutional layers up to the layer in the network under study. A weighted combination of these feature maps gives the *importance score* $M_c^{(i,j)} \in [0, 1]$ that indicates the importance of the feature in the image at location $(i, j)$ for the class cth. This score is given as

$$M_c^{(i,j)} = \frac{1}{Z} \sum_k w_k \frac{\partial y_c}{\partial F_k^{(i,j)}}, \qquad (B1)$$

where $Z = \sum_k w_k$ is the normalization factor, corresponding to the sum of the $k$th weights, while $w_k$ is the weight corresponding to the feature map $F_k$. In our case, we focus on the neutral regions, categorized with a value of $c = 1$ in our maps. Values close to one indicate high importance, while in the opposite case, it indicates irrelevance.

In Fig. B1, we show the result for three hidden layers. The first column shows the input image and the $M_c$ score represented by dark shadows, indicating the location in the image employed in the classification of the neutral regions. Solid line contours indicate the ground truth. The central column shows the Grad-CAM filtered

region obtained by element-wise multiplication of the input image with the importance score. This plot shows us what features the network emphasizes in the image for identifying the neutral region. The right column visualizes the hidden layer output, with the number of subpanels corresponding to the number of channels. From top to bottom, we have the output of the convolution block at the second level of the encoder after two convolutional layers and a pooling operation. The hidden state has angular and channel dimensions (64, 64, 32). We can see that in the encoder, the network focuses on the regions with the highest intensity, which, thanks to the pre-process presented in Section 3, are mostly located within the neutral region. The different channels in the hidden layer show a similar conclusion, with the convolutional operation capturing the large-scale region that produces 21-cm signals. In the second row, the bottom of the U-Net, known as the low-dimensional latent space, with dimension (16,16,128), gives a compressed representation of the input image, and it appears to focus on location in the image with the highest and lowest values. Our interpretation is that the network focuses on these extreme values to quantify the ‘threshold’ value that sets the boundary between neutral and ionized regions. This interpretation is also supported by the 128 hidden layer plots in the right-hand panel, as the compressed data shows different constant values across the channels. In the last row, we show the importance score from the final convolution before the binarization of the output. Here, it appears that the network uses the threshold value defined in the bottom layer of the U-Net to locate the delimitation that defines the neutral regions, as the shadow is located along the contour of the ground truth (black solid line). We notice that some locations in the image with substantial foreground residuals are wrongly included. The hidden layer plot shows that the network struggles to remove the foreground residual completely.
<table><caption>**Table A1.** SegU-Net hyperparameter optimization analysis for the best-performing set-ups of seven parameters with optimization on the validation loss.</caption><thead><tr><th>Ranking</th><th>Activation</th><th>Channels latent space</th><th>Depth</th><th>Dropout</th><th>Final activation</th><th>Kernel</th><th>Pooling type</th><th>r<sub>φ</sub> (per cent)</th><th>Validation loss (× 10<sup>-2</sup>)</th></tr></thead><tbody><tr><td>1</td><td>LeakyReLu</td><td>256</td><td>3</td><td>0.42</td><td>σ(x)</td><td>6</td><td>max</td><td>89.08</td><td>6.59</td></tr><tr><td>2</td><td>LeakyReLu</td><td>128</td><td>4</td><td>0.00</td><td>σ(x)</td><td>5</td><td>max</td><td>89.02</td><td>6.62</td></tr><tr><td>3</td><td>Elu</td><td>128</td><td>3</td><td>0.34</td><td>x</td><td>11</td><td>average</td><td>87.76</td><td>7.27</td></tr><tr><td>4</td><td>LeakyReLu</td><td>128</td><td>4</td><td>0.50</td><td>σ(x)</td><td>5</td><td>max</td><td>88.72</td><td>6.85</td></tr><tr><td><b>5</b></td><td><b>ReLu</b></td><td><b>256</b></td><td><b>4</b></td><td><b>0.05</b></td><td><b>x</b></td><td><b>7</b></td><td><b>average</b></td><td><b>88.50</b></td><td><b>7.48</b></td></tr><tr><td>6</td><td>ReLu</td><td>256</td><td>5</td><td>0.14</td><td>x</td><td>7</td><td>max</td><td>86.53</td><td>9.15</td></tr></tbody></table>

## APPENDIX B: INSIDE THE BLACK BOX

The trained model we presented in this paper is able to recover the ionized field from noisy images with residual foreground contamination. This is an indication that the network learns to identify the regions of interest from important hidden features that maximize the recovery. However, the machine learning model's complexity, high dimensionality and non-linearity make them difficult to interpret and regulate, so these applications are often referred to as a *black box*. Here, we present a first attempt to open and look inside SegU-Net black box. A standard tool to visualize and understand the decisions made by a general convolutional neural network is the Gradient-weighted Class Activation Mapping (Grad-CAM) technique (Selvaraju et al. 2019). This method applied in segmentation and object classification highlights the features of an input image employed in predicting a particular class. Grad-CAM achieves this by computing the gradients of the predicted class $y_c \equiv x_{\text{HI}}^B$ score with respect to the feature maps $F_k$ of all the $k > 0$ convolutional layers up to the layer in the network under study. A weighted combination of these feature maps gives the *importance score* $M_c^{(i,j)} \in [0, 1]$ that indicates the importance of the feature in the image at location $(i, j)$ for the class cth. This score is given as

$$M_c^{(i,j)} = \frac{1}{Z} \sum_k w_k \frac{\partial y_c}{\partial F_k^{(i,j)}}, \qquad (B1)$$

where $Z = \sum_k w_k$ is the normalization factor, corresponding to the sum of the $k$th weights, while $w_k$ is the weight corresponding to the feature map $F_k$. In our case, we focus on the neutral regions, categorized with a value of $c = 1$ in our maps. Values close to one indicate high importance, while in the opposite case, it indicates irrelevance.

In Fig. B1, we show the result for three hidden layers. The first column shows the input image and the $M_c$ score represented by dark shadows, indicating the location in the image employed in the classification of the neutral regions. Solid line contours indicate the ground truth. The central column shows the Grad-CAM filtered

region obtained by element-wise multiplication of the input image with the importance score. This plot shows us what features the network emphasizes in the image for identifying the neutral region. The right column visualizes the hidden layer output, with the number of subpanels corresponding to the number of channels. From top to bottom, we have the output of the convolution block at the second level of the encoder after two convolutional layers and a pooling operation. The hidden state has angular and channel dimensions (64, 64, 32). We can see that in the encoder, the network focuses on the regions with the highest intensity, which, thanks to the pre-process presented in Section 3, are mostly located within the neutral region. The different channels in the hidden layer show a similar conclusion, with the convolutional operation capturing the large-scale region that produces 21-cm signals. In the second row, the bottom of the U-Net, known as the low-dimensional latent space, with dimension (16,16,128), gives a compressed representation of the input image, and it appears to focus on location in the image with the highest and lowest values. Our interpretation is that the network focuses on these extreme values to quantify the 'threshold' value that sets the boundary between neutral and ionized regions. This interpretation is also supported by the 128 hidden layer plots in the right-hand panel, as the compressed data shows different constant values across the channels. In the last row, we show the importance score from the final convolution before the binarization of the output. Here, it appears that the network uses the threshold value defined in the bottom layer of the U-Net to locate the delimitation that defines the neutral regions, as the shadow is located along the contour of the ground truth (black solid line). We notice that some locations in the image with substantial foreground residuals are wrongly included. The hidden layer plot shows that the network struggles to remove the foreground residual completely.
<figure><img src="image_11.png" /><figcaption>**Figure B1.** Region of interest detected by Grad-CAM for three hidden layers. *Left-hand panels:* Input image with shadow areas that indicate the region of attention detected by the Grad-CAM method. Black solid contours indicate the ground truth for comparison. *Central panel:* The filtered Grad-CAM image element-wise multiplication between the input and the $M_c$ filter. *Right-hand panel:* A visualization of the hidden layer output. The number of subpanels indicates the channel size of the hidden layer.</figcaption></figure>

This paper has been typeset from a TeX/LATEX file prepared by the author.
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Region of interest detected by Grad-CAM for three hidden layers. *Left-hand panels:* Input image with shadow areas that indicate the region of attention detected by the Grad-CAM method. Black solid contours indicate the ground truth for comparison. *Central panel:* The filtered Grad-CAM image element-wise multiplication between the input and the $M_c$ filter. *Right-hand panel:* A visualization of the hidden layer output. The number of subpanels indicates the channel size of the hidden layer.</figcaption></figure>

This paper has been typeset from a TeX/LA TeX file prepared by the author.