#  Bayesian Data Fusion of Imperfect Fission Yields for Augmented Evaluations

Z.A. Wang  State Key Laboratory of Nuclear Physics and Technology, School of
Physics, Peking University, Beijing 100871, China  J.C. Pei<EMAIL_ADDRESS>State Key Laboratory of Nuclear Physics and Technology, School of Physics,
Peking University, Beijing 100871, China  Y.J. Chen  China Nuclear Data
Center, China Institute of Atomic Energy, Beijing 102413, China  C.Y. Qiao
State Key Laboratory of Nuclear Physics and Technology, School of Physics,
Peking University, Beijing 100871, China  F.R. Xu  State Key Laboratory of
Nuclear Physics and Technology, School of Physics, Peking University, Beijing
100871, China  Z.G. Ge  China Nuclear Data Center, China Institute of Atomic
Energy, Beijing 102413, China  N.C. Shu  China Nuclear Data Center, China
Institute of Atomic Energy, Beijing 102413, China

###### Abstract

We demonstrate that Bayesian machine learning can be used to treat the vast
amount of experimental fission data which are noisy, incomplete, discrepant,
and correlated. As an example, the two-dimensional cumulative fission yields
(CFY) of neutron-induced fission of 238U are evaluated with energy
dependencies and uncertainty qualifications. For independent fission yields
(IFY) with very few experimental data, the heterogeneous data fusion of CFY
and IFY is employed to interpolate the energy dependence. This work shows that
Bayesian data fusion can facilitate the further utilization of imperfect raw
nuclear data.

Accurate modeling of nuclear fission remains an open problem although its
significant applications are generally known future . In fact, for modern
nuclear energy productions in a more compact and sustainable, cleaner and
safer way, a supply of high quality data of multiple fission observables as a
function of incident neutron energy is crucial nd . In addition, more accurate
fission yields are needed to address the reactor anti-neutrino anomaly RAA ;
RAA2 . In major nuclear data libraries endf ; jendl ; jeff ; cendl , evaluated
fission yields are only available for neutron incident energy at thermal
energy, 0.5 and 14 MeV. Both experimental measurements and reliable modelings
of nuclear fission are notoriously difficult. Nevertheless, there was a vast
amount of raw experimental fission data exfor  but these data are generally
noisy, incomplete, and discrepant. Therefore it is desirable to exploit the
further value of these imperfect fission data for novel nuclear applications
as well as for constraints on fission theories.

The spreading fission product yields (FPY) are the most frequently studied
observables, and FPY distributions evolve in time due to cascade decays of
fragments. There are more data of cumulative fission yields (CFY) which are
measured after post-fission $\beta$-decays but much fewer data of independent
fission yields (IFY) which are identified promptly after fission exfor .
Conventionally, semi-empirical models such as Brosa brosa  and GEF gef  models
with physics guidances are used for evaluation of fission data. However, semi-
empirical models invoke many parameters, which undermines their
applicabilities for new data and uncertainty quantifications. Advanced
evaluation methods should take into account the correlations between different
fission data cgmf . In this context, the covariance matrix is widely used for
uncertainty assessments endf ; cov2 ; cov3 .

Machine learning is now revolutionizing all scientific fields and is
particularly powerful for treatments of complex big data. The belief behind
data science is that data includes all correlations. There have been
increasing applications of machine learning in nuclear physics (see a review
ai ). In particular, machine learning can be used to extract and discover
knowledge from existing nuclear data, which is referred to as data mining. For
example, machine learning has been applied for inferences of nuclear
structures wn-m ; utama-1 ; zmniu ; bai ; utama ; xu  and nuclear reactions ma
; alovell21 ; Neudecker ; wang . The application of machine learning to
augment evaluations of nuclear data is very anticipated augmented . Previously
we have applied the Bayesian neural network (BNN) for evaluation of incomplete
fission mass yields wangza1  and discrepant charge yields qiao . The multi-
layer BNN has been optimized for evaluation of fission yields with physical
odd-even effects qiao  and non-negative constraints wangza2 . There was also
machine learning of fission yields using the mixture density network method
alovell20 . Besides, the Gaussian process can also solve regression problems
of nuclear data evaluations but is assumed with local correlations.

In this Letter, we aim to make the best utilization of vast raw experimental
fission data by data fusion fusion , which is a prevalent way to deal with
imperfect data and includes all underlying correlations. The fused data is
expected to be more informative than separated data sources. It is possible
that some non-local and high-dimensional correlations could be weak but non-
negligible, in analogy to long-range and many-body interactions. The inference
would be less precise when data in some energies is sparse, however, its
correlations with other data in other energies can be beneficial to improve
the inference. It is known that Bayesian machine learning is ideal for
uncertainty quantification brand  which is key in fusion of imperfect data.
Furthermore, the fusion of complex heterogeneous data can be naturally treated
by Bayesian machine learning, which exhibits the powerful capabilities of
machine learning by including correlations beyond the Gaussian process and the
covariance matrix.

We employ BNN bnn  for machine learning of distributions of fission yields.
The posterior distribution $p(\omega|{\rm x},t)$ of BNN is based on a prior
distribution $p(\omega)$ of network parameters $\omega$ and a likelihood
function $p({\rm x},t|\omega)$,

$p(\omega|{\rm x},t)=\frac{p({\rm x},t|\omega)p(\omega)}{p({\rm x},t)}.$ (1)

Thus the resulting inference has a distribution and naturally provides the
associated uncertainty. The data set is given as $D=\\{x_{i},t_{i}\\}$, where
$x_{i}$ is the input and $t_{i}$ is the output fission yields. The details of
BNN has been described in previous works wangza1 ; brand . For the data
fusion, it is crucial to taken into account the experimental uncertainties by
the likelihood function,

$p({\rm x},t\mid\omega)=\exp(-\chi^{2}/2).$ (2)

The cost function $\chi$2($\omega$) reads:

$\chi^{2}(\omega)=\sum^{N}_{i=1}\frac{(t_{i}-f({\rm
x_{i}},\omega))^{2}}{\delta_{i}^{2}+\sigma_{i,\rm expt}^{2}},$ (3)

where $f({\rm x_{i}},\omega)$ denotes the network values. Note that the
weights in the likelihood include a hyperparameter noise scale
$\delta_{i}^{2}$ and the experimental uncertainty $\sigma_{i,\rm expt}^{2}$.
The hyperparameter $\delta_{i}^{2}$ is uniform for all data points and changes
in the learning process until numerical convergence. It is reasonable to see
that data points with large experimental uncertainties would have small
weights in the data fusion. It is also interesting to study how the
experimental uncertainties prorogate to the prediction uncertainties. The data
fusion is in some sense similar to the conception of model averaging and model
mixing wn-m  but doesn’t suffer from the bias in model selections.

Firstly, we demonstrate the uncertainty prorogations by varying experimental
uncertainties, or by deleting or adding data points deliberately to mimic
incomplete and discrepant data. We use BNN to evaluate the fission mass
distribution of 238U as an illustrative example, with experimental data from
Ref. Ramos . The data set also includes the evaluated independent fission
yields of 238U with neutron incident energies at 0.5 and 14 MeV from JENDL
jendl . Note that neutron-induced fission of 238U doesn’t occur at 0.5 MeV and
the yield data of 0.5 MeV is a hypothetical evaluation. Here BNN adopts a
single hidden layer with 10 neurons because the tests employ a small data set
of 309 data points. Fig.1(a) displays the evaluation of fission yields of 238U
by including original experimental uncertainties. The BNN evaluations
generally agree with experimental data except that the evaluation at mass
$A$=100 is lower than experiments, since the JENDL evaluation is also lower.
Fig.1(b) shows the inferred uncertainties by varying the experimental
uncertainties in Fig.1(a) with a factor from 0.5 to 2.0. We see that
uncertainties indeed increases with increasing experimental uncertainties but
don’t increase linearly. It is reasonable to see that the uncertainties are
generally larger around two peaks. Fig.1(c) and Fig.1(d) show the predictions
when the experimental data are deleted at the right or left peak,
respectively. We clearly see that the resulted uncertainty of the absent data
is larger than that with experimental data. Fig.1(e) and Fig.1(f) show the
influences of discrepant data. In Fig.1(e), we add some extra data at the left
peak that are not far from experiments with small uncertainties. We see the
the evaluation moves toward the extra data with increased uncertainties. In
Fig.1(f), we add some extra data that are far from experiments but with large
uncertainties, However, we see the evaluation doesn’t move significantly
towards the extra data. Therefore the uncertainty prorogation is a
comprehensive effect and would not significantly change due to a specific
data.

Figure 1:  The BNN evaluaitons of fission yields from neutron-induced fission
of 238U, with experimental data from Ramos . The shadows denote the BNN
uncertainties as given by the confidential interval (CI) at 95$\%$. (a) shows
the BNN evaluation with experimental uncertainties. (b) shows the evaluated
uncertainties by varying experimental uncertainties with a factor from 0.5 to
2.0. (c) and (d) shows the evaluations when data points of the right peak or
the left peak are deleted, respectively. (f) and (g) show the evaluations with
some extra noisy data (red color).

Next we did the practical Bayesian data fusion of CFY of neutron-induced
fission of 238U from different experiments with different incident energies.
The raw imperfect experimental fission yields are taken from the EXFOR library
exfor ; exfor1 . In BNN, input variables are given in terms of ($Z$, $N$,
$E$), i.e., the atomic number $Z$, neutron number $N$ of fragments and the
neutron incident energy $E$. There are total 33 experiments with different
energies of CFY of 238U with about 1221 scattering data points. In addition,
2064 data points of evaluated CFY from JENDL at energies of 0.5 and 14 MeV are
included as a learning constraint. Here BNN adopts a double-layer network with
20-20 neurons.

Figure 2:  The BNN evaluated yield-energy relations for 99Mo (a), 135Xe (b),
140Ba (c) and 147Nd (d) from $n$+238U fission are shown. The raw experimental
data are taken from EXFOR exfor . Different experiments are denoted by
different colours. The shadows denote the corresponding BNN uncertainties
given by CI at 95$\%$.  Figure 3:  The two-dimensional CFY distributions of
$n$+238U fission obtained by BNN data fusion. Panels (a), (b) and (c) shows
CFY at incident neutron energies of 2, 8 and 14 MeV respectively. Panels (d),
(e) and (f) shows the corresponding uncertainties at different energies.

The yield-energy relations of some long-lived isotopes are of particular
application interests for monitoring fission environments. For example, 135Xe
has large neutron absorption cross sections and is a “poison” for reactors 135
. Fig.2 shows the evaluated energy dependence of fission yields of 99Mo,
135Xe, 140Ba and 147Nd. In Fig.2, we see that the experimental data are indeed
noisy, sparse and even discrepant. We see that data fusion can give reasonably
the yield-energy relations and uncertainty quantifications. Since JENDL
evaluations are included, consequently BNN results are close to JENDL jendl
and ENDF endf  evaluations at 0.5 and 14 MeV for the above 4 fragment yields.
Our key motivation is to infer the energy dependence between 0.5 and 14 MeV.
For 99Mo, CFY decreases firstly but doesn’t decrease after 10 MeV. For 140Ba,
CFY decreases smoothly as energy increases. There are more data for 99Mo and
140Ba and they are usually adopted as standards to determine other yields. For
135Xe, there are very few experimental data and resulting uncertainties are
larger than other cases. For 147Nd, there are some significantly discrepant
data between 5$\sim$10 MeV. Our BNN evaluations are flat before 10 MeV and are
close to the latest measurement gooden . The corresponding uncertainties
become larger in this energy regime. In Fig.2, the inference uncertainties are
larger than some experimental uncertainties but they are acceptable regarding
the influences of discrepant data. The extrapolated yields at higher energies
generally have increasing uncertainties. The BNN uncertainty includes two
parts: the overall regression noise and the data-dependent posterior
uncertainty bnn , while they are not exclusive. In Fig.2, there is an overall
noise scale of $\sim$0.6 which is related to the description capability of BNN
and can be reduced by more complicated neural network or physics-informed
neural network. Besides, Fig.1 and Fig.2 demonstrate that the experimental
uncertainty, data sparsity and data discrepancy are reflected by the data-
dependent uncertainty. Usually different experiments adopt different
criterions for normalization of fission yields. The BNN evaluations can be
further improved by careful calibrations of raw experimental data.

Fig.3 displays the two-dimensional CFY distributions of neutron-induced
fission of 238U at energies of 2, 8 and 14 MeV from BNN data fusion. Recent
experiments have made great progress in measurements of the isotopic mass
yields using the inverse kinetic method Ramos , and it is still very difficult
to obtain the complete distributions. The two-dimensional CFY distribution is
very different from Gaussian distributions and it is difficult to be evaluated
directly by semi-empirical models. In principle, our results include the
yield-energy relations of all fragments. We can see that the highest yields
are around $Z$=53$\thicksim$54 and $N$=79$\thicksim$81, which decrease
gradually with increasing energies. It is known that symmetric fission would
become more prominent as excitation energy increases symf1 ; symf2 . The two
peaks associated with asymmetric fission modes reduce with increasing
energies, while the central distributions associated with the symmetric
fission mode are increasing. Therefore, the energy dependence of fission modes
can be reasonably described by BNN data fusion. Fig.3 (d, e, f) shows the
corresponding uncertainties. There is an overall background noise scale
$\sim$0.6 for all data. Obviously the uncertainty is energy dependent. The
average uncertainties of Fig.3 (d, e, f) are 0.789, 0.805 and 0.793 for 1032
points, respectively. The larger uncertainty at 8 MeV is related to the data
sparsity around this energy regime. The uncertainties of some FPY are as large
as 1.4 where CFY have large values. The relative uncertainties around peaks
are actually smaller. The one-dimensional charge yields can also be extracted
from the two-dimensional CFY distributions. The precise and complete CFY are
crucial to estimate the abnormal anti-neutrino spectrum at reactors RAA2  and
our results will be useful for such studies.

Figure 4:  The two-dimensional IFY distribution of $n$+238U fission with
neutron incident energy at 8 MeV obtained by BNN data fusion. (a) shows the
inference by learning pure IFY data. (b) shows the inference by heterogeneous
data fusion of CFY and IFY. (c) and (d) shows the uncertainties corresponding
to (a) and (b), respectively.

Different from CFY, there are much fewer experimental data of IFY and its
evaluation is more challenging. IFY distributions are particularly useful for
constraints of fission theories. For 238U, there are only 211 IFY experimental
data exfor  and 2064 evaluated data points from JENDL. For better
interpolation of energy dependence of IFY, we employ the heterogeneous data
fusion of CFY and IFY, by utilizing the correlations between CFY and IFY. The
learning dataset is ($A_{i}$, $Z^{C/I}_{i}$, $E_{i}$, Y${}^{C/I}_{i}$), in
which CFY and IFY share the mass number $A$ of fragments and the incident
energy $E$, but their proton numbers are deliberately separated into two
groups as $Z^{I}_{i}$ and $Z^{C}_{i}$. Such a scheme is based on the fact that
CFY mainly results from $\beta$-decays of IFY but they have very different
distributions. The heterogeneous datasets are learned by the same network so
that their multi-dimensional correlations are naturally included. The
evaluated IFY from JENDL at 0.5 MeV and 14 MeV are used in the learning and
our main purpose is to interpolate the energy dependence at other energies.
Fig.4 displays the evaluation of IFY at 8 MeV with learning of pure IFY and
heterogeneous CFY and IFY, respectively. The homogeneous fusion of IFY adopts
a 16-16 network and the heterogeneous fusion of CFY and IFY adopts a 20-20
network. There are almost no experimental IFY data between 2 and 14 MeV. It is
known that the energy dependence of PFY is nonlinear from Fig.2. Thus the
evaluation of IFY at 8 MeV is not reliable by learning only IFY data. We see
that there are almost no FPY from the symmetric fission between two peaks in
Fig. 4(a), but there are considerable symmetric fission FPY in Fig. 4(b). This
demonstrated that the heterogeneous data fusion of CFY and IFY can build the
energy dependence information of CFY into IFY. The resulting IFY in Fig. 4(b)
are consistent with the CFY in Fig.3(b). Fig. 4(c)(d) displays the
uncertainties correspondingly. We see that the homogeneous fusion has a much
smaller overall noise but some points have large uncertainties. The
heterogeneous fusion has an overall background noise scale about 0.2 due to
the influence of CFY. Indeed, it has been pointed out that cross-experiment
correlations would result in increased final estimated uncertainties cov3 . On
the other hand, without the background noise, the data-dependent uncertainties
of heterogeneous data fusion are generally smaller than that of the
homogeneous fusion.

In summary, we have applied Bayesian machine learning for the evaluation of
imperfect fission yields, in which the raw experimental data are generally
noisy, incomplete, discrepant and correlated. The Bayesian data fusion can be
used to augment the inference of energy dependence of imperfect fission data
by including underlying correlations, which is a crucial need for next-
generation nuclear energy productions. The two-dimensional distributions of
cumulative fission yields in terms of energy dependence are obtained. The
evolution of fission modes is revealed by the smooth transition of two-
dimensional fission yields. The yield-energy relations of key fragments are
now given with uncertainty quantifications, in which energy dependencies of
uncertainties are data dependent. We also applied the heterogeneous data
fusion to interpolate the energy dependence of independent fission yields
which have very few experimental data. We expect this method can also be
applied to the heterogeneous data fusion of multiple fission observables and
flexible physics constraints for a comprehensive understanding of nuclear
fission.

###### Acknowledgements.

We are grateful for valuable comments by W. Nazarewicz. This work was
supported by the National Key R$\&$D Program of China (Contract No.
2018YFA0404403) and by the National Natural Science Foundation of China under
Grants No. 11975032, 11835001, 11790325, and 11961141003.

## References

  * (1) M. Bender et al., J. Phys. G 47, 113002(2020).
  * (2) L.A. Bernstein, D. A. Brown, A. J. Koning, B.T. Rearden, C. E. Romano, A. A. Sonzogni, A. S. Voyles, and W. Younes, Ann. Rev. Nucl. Part. Sci. 69, 109(2019).
  * (3) G. Mention, M. Fechner, Th. Lasserre, Th. A. Mueller, D. Lhuillier, M. Cribier, and A. Letourneau, Phys. Rev. D 83, 073006 (2011).
  * (4) A.A. Sonzogni, E.A. McCutchan, T.D. Johnson, and P. Dimitriou, Phys. Rev. Lett. 116, 132502(2016).
  * (5) M.B. Chadwick, et al., Nuclear Data Sheets 112, 2887 (2011).
  * (6) K. Shibata, et al., J. Nucl. Sci. Tech. 48, 1(2011).
  * (7) Joint Evaluated Fission and Fusion (JEFF) Nuclear Data Library, https://www.oecd-nea.org/dbdata/jeff/
  * (8) Z.G. Ge, Z. X. Zhao, H. H. Xia, Y. X. Zhuang, T. J. Liu, J. S. Zhang and H. C. Wu, J. Korean Phys. Soc. 59, 1052 (2011)
  * (9) https://www-nds.iaea.org/exfor/
  * (10) U. Brosa, S. Grossmann, A. Müller, Phys. Rep. 197, 167 (1990).
  * (11) K.H. Schmidt, B. Jurado, C. Amouroux, and C. Schmitt, Nuclear Data Sheets 131, 107(2016).
  * (12) P. Talou, R. Vogt, J. Randrup, et al., Eur. Phys. J. A 54, 9 (2018).
  * (13) K. Tsubakihara, S. Okumura, C. Ishizuka, T. Yoshida, F. Minato, and S. Chiba, J. Nucl. Sci. Tech. 58, 151(2021).
  * (14) P. Talou, EPJ Nuclear Sci. Technol. 4, 29 (2018).
  * (15) P. Bedaque, A. Boehnlein, M. Cromaz, et al., Eur. Phys. J. A 57, 100 (2021)
  * (16) L. Neufcourt, Y. Cao, S. A. Giuliani, W. Nazarewicz, E. Olsen, and O. B. Tarasov, Phys. Rev. C 101, 044307 (2020).
  * (17) R. Utama, J. Piekarewicz, and H. B. Prosper, Phys. Rev. C 93, 014311 (2016).
  * (18) Z. M. Niu and H. Z. Liang, Phys. Lett. B 778, 48 (2018).
  * (19) R. Utama, W.C. Chen, and J. Piekarewicz, J. Phys. G 43, 114002 (2016).
  * (20) D. Wu, C. L. Bai, H. Sagawa, and H. Q. Zhang, Phys. Rev. C 102, 054323 (2020).
  * (21) Y. Ma, C. Su, J. Liu, Z. Ren, C. Xu, and Y. Gao, Phys. Rev. C 101, 014304 (2020).
  * (22) C. W. Ma, D. Peng, H. L. Wei, Z. M. Niu, Y. T. Wang, and R. Wada, Chin. Phys.C 44, 014104 (2020).
  * (23) A. E. Lovell, F. M. Nunes, M. Catacora-Rios, and G. B. King, J. Phys. G 48, 014001 (2021).
  * (24) D.Neudecker, O.Cabellos, A. R. Clark, M. J. Grosskopf, W. Haeck, M. W. Herman, J. Hutchinson, T. Kawano, A. E. Lovell, I. Stetcu, P. Talou, and S. V. Wiel, Phys. Rev. C 104, 034611(2021).
  * (25) Y. Wang, F. Li, Q. Li, H. Lu, K. Zhou, Phys. Lett. B 822,136669(2021).
  * (26) P. Vicente Valdez, L. Bernstein, M. Fratoni, Ann. Nucl. Energy 163, 108596(2021).
  * (27) Z. A. Wang, J. C. Pei, and Y. Liu, Y. Qiang, Phys. Rev. Lett. 123, 122501 (2019).
  * (28) C. Y. Qiao, J. C. Pei, Z. A. Wang, Y. Qiang, Y. J. Chen, N. C. Shu, and Z. G. Ge, Phys. Rev. C 103, 034621(2021).
  * (29) Z.A. Wang, J.C. Pei, arXiv:2106.11746, Phys. Rev. C (accepted).
  * (30) A. E. Lovell, A. T. Mohan, and P. Talou, J. Phys. G 47, 114001 (2020).
  * (31) T. Meng, X. Jing, Z. Yan, W. Pedrycz, Information Fusion 57, 115(2020).
  * (32) D.R. Phillips, R.J. Furnstahl, U. Heinz, T. Maiti, W. Nazarewicz, F.M. Nunes, M. Plumlee, M.T. Pratola, S. Pratt, F.G. Viens, and S.M. Wild, J. Phys. G 48, 072001(2021).
  * (33) R. Neal, Bayesian Learning of Neural Network, Springer, New York (1996).
  * (34) D. Ramos et al., Phys. Rev. C 101, 034609(2020).
  * (35) N.Otuka et al., Nucl. Data Sheets 120, 272(2014).
  * (36) O. Moreira, Ann. Nucl. Energy 39, 62(2012).
  * (37) M.E. Gooden, et al., Nuclear Data Sheets 131, 319(2016)
  * (38) J. A. Sheikh, W. Nazarewicz, and J. C. Pei, Phys. Rev. C 80, 011302(R)(2009).
  * (39) J. Zhao, T. Nikšić, D. Vretenar, S.G.Zhou, Phys. Rev. C 99, 014618 (2019).


# ReservoirComputing.jl: An Efficient and Modular Library for Reservoir
Computing Models

Francesco Martinuzzi
Remote Sensing Centre for Earth System Research
Center for Scalable Data Analytics and Artificial Intelligence
Julia Computing Inc.
<EMAIL_ADDRESS>
Chris Rackauckas
Massachusetts Institute of Technology
Julia Computing Inc.
<EMAIL_ADDRESS>
Anas Abdelrehim
Julia Computing Inc.
<EMAIL_ADDRESS>
Miguel D. Mahecha
Remote Sensing Centre for Earth System Research
<EMAIL_ADDRESS>
Karin Mora
Remote Sensing Centre for Earth System Research
German Centre for Integrative Biodiversity Research (iDiv)
<EMAIL_ADDRESS>

###### Abstract

We introduce ReservoirComputing.jl, an open source Julia library for reservoir
computing models. The software offers a great number of algorithms presented
in the literature, and allows to expand on them with both internal and
external tools in a simple way. The implementation is highly modular, fast and
comes with a comprehensive documentation, which includes reproduced
experiments from literature. The code and documentation are hosted on Github
under an MIT license https://github.com/SciML/ReservoirComputing.jl.

_Keywords_ reservoir computing, echo state networks, julia, machine learning

##  1 Introduction

Time series modeling is a very common technique throughout many areas of
machine learning. However, many standard recurrent models are known to be
susceptible to problems such as the vanishing gradient (Pascanu et al., 2013)
or the extreme sensitivity of chaotic systems to their parameterization
(Wiggins et al., 2003). To counter these issues reservoir computing (RC)
techniques were introduced as recurrent models which can be trained without
requiring gradient-based approaches (Lukoševičius and Jaeger, 2009).
Independently proposed as echo state networks (ESNs) (Jaeger, 2001) and liquid
state machines (LSMs) (Maass et al., 2002), these architectures are based on
the expansion of the input data using a fixed random internal layer, known as
the reservoir, and the subsequent mapping of the reservoir to match an output.
This output mapping can normally be written in a simple closed form, such as
an $L_{2}$ minimization computed via a single QR-factorization. This allows
the RC approach to achieve faster computational times with less parameter
tuning compared to deep learning models which rely on local optimization.
Numerous alterations have been made since their inception, such as deep
architectures (Gallicchio et al., 2018) and different training approaches
Chatzis and Demiris (2011); Shi and Han (2007). Applications have confirmed
that these improved architectures and training techniques can be impactful in
many fields of study: from surrogates for stiff systems (Anantharaman et al.,
2020) to reconstruction of chaotic attractors (Pathak et al., 2018a) and
prediction of climate dynamics (Nadiga, 2021). This expanding landscape
necessitates the development of production-quality software to empower
scientists with all of the latest techniques without having to implement from
scratch every novel variation on their own. While a number of libraries for
ESNs and RC are available (Trouvain et al., 2020; Steiner et al., 2021;
Pontes-Filho et al., 2020), they lack the hackability necessary to alter the
provided architectures and achieve the top performances from the literature.
To provide a truly modular library for RC models we present
ReservoirComputing.jl, an efficient and flexible library written using the
Julia language (Bezanson et al., 2017). This software provides a user oriented
package with intuitive high level application programming interfaces (APIs),
while maintaining a great level of customization. This flexibility is obtained
from both the design choices of ReservoirComputing.jl and by leveraging the
multiple dispatch architecture of Julia. This gives users the freedom to mix
any Julia code with the reservoir computers in a seamless way to promote the
latest performance and research.

##  2 Theoretical Background

The setup of a RC model (Konkoli, 2018) is based on leveraging the complex
dynamics of a nonlinear system, the reservoir $\mathcal{R}$. The role of this
layer is to expand the input data $\mathbf{u}(t)$ at time $t$ into a higher
dimension. This operation can be aided by an input layer $\chi$, which maps
the data onto the reservoir. The input data $\mathbf{u}(t)$ drives the
reservoir $\mathcal{R}$ dynamics to a certain configuration. The resulting
state $\mathbf{x}(t)$ represents the expansion of the input data and is used
for training. Once the states have been obtained for all the training length
$T$ the RC model can be trained. The states $\mathbf{x}(t)$ can be manipulated
at this stage: through concatenation with the input data $\mathbf{u}(t)$ such
as $\mathbf{z}(t)=[\textbf{x}(t),\textbf{u}(t)]$ (Jaeger, 2007), using padding
with some constant value $c$ usually set $c=1.0$ such as
$\mathbf{z}(t)=[c,\textbf{u}(t)]$ (Lukoševičius, 2012), or through nonlinear
transformations (Chattopadhyay et al., 2020). The training is performed in one
step in order to find the best fit between the states $\mathbf{x}(t)$ and the
desired output $\mathbf{y}(t)$g. The most common way to perform the training
is linear regression, creating the readout layer $\psi$. This layer is finally
used in the prediction phase to obtain the desired output
$\mathbf{v}(t)=\psi(\mathbf{x}(t))$.

##  3 Library Overview

Figure 1: ReservoirComputing.jl architecture. _SV_ stands for support vector.

ReservoirComputing.jl features the implementation of different RC models,
training methods, and state variations. The architecture with the main
components are detailed below and illustrated in Figure 1. Given the modular
nature of the software more algorithms can be easily implemented on top of the
existing ones. Any of the many libraries of the Julia community or custom
codes can be used as part of the functions within the reservoir computer.
Importantly, as Julia is a just-in-time (JIT) compiled language, these
extensions do not sacrifice performance.

### Building models.

ReservoirComputing.jl provides a simple interface for model building that
closely follows the workflow presented in the corresponding literature. It
includes a standard implementation of ESNs (Lukoševičius, 2012) as well as a
hybrid variation (Pathak et al., 2018b), gated recurrent unit ESN (Wang et
al., 2020) and double activation function ESN (Lun et al., 2015). Multiple
input layers $\chi$ and reservoirs $\mathcal{R}$ are also provided, ranging
from weighted input layers (Lu et al., 2017) to minimally complex input layers
and reservoirs (Rodan and Tino, 2010; Rodan and Tiňo, 2012), including a
reservoir obtained through pseudo single value decomposition (Yang et al.,
2018). Reservoir computing with cellular automata (ReCA) (Yilmaz, 2014;
Nichele and Molund, 2017) is another family of models available in the library
leveraging the package CellularAutomata.jl (Martinuzzi, 2022).

### Training.

The training algorithms are implemented in a low level fashion, supporting all
RC models of the library. Multiple training methods to obtain the output layer
$\psi$ can be obtained from open source libraries such as MLJLinearModels.jl
(Blaom and Vollmer, 2020), GaussianProcesses.jl (Fairbrother et al., 2021) and
LIBSVM.jl (Kornblith and Pastell, 2021), a Julia porting of LIBSVM (Chang and
Lin, 2011).

### Prediction.

ReservoirComputing.jl leverages diverse prediction techniques. The generative
approach allows the model to run autonomously: the predicted output
$\mathbf{v}(t)$ is fed back into the model to obtain the prediction for the
next time step $\mathbf{v}(t+1)$. The predictive approach instead uses the
standard feature-label setup.

### States modifiers.

Other lower level implementations are included, such as the possibility to
modify the state vectors. All the approaches detailed in §2 are included in
the library.

##  4 Code Quality

The library is made available on Github, where it is continuously tested
through Github actions. The package is part of the scientific machine learning
(SciML) community (sci, ), a NumFOCUS sponsored project since 2020 (num, ). An
extensive online documentation is also provided with the library, where models
are documented at length with multiple examples for different applications.

A comparison of features and performance of ReservoirComputing.jl with similar
libraries is provided in Table 1 and in Figure 2, respectively. The feature
comparison in Table 1 showcases the large offering of models available. The
plot in Figure 2 illustrates the better computational speed of the library.
The performance test task is a next step prediction of the Mackey-Glass system
(Glass and Mackey, 2010) with time delay $\tau=17$. The dense reservoir matrix
and the dense input matrix are generated with uniform distribution sampled
from $[-1,1]$. The spectral radius of the reservoir matrix is scaled by 1.25.
The ridge regression parameter is set to $10^{-8}$. Training and prediction
lengths are both equal to 4999. The time reported is the sum of training and
prediction. For central processing unit (CPU) computations the precision is
set to float64, for the graphics processing unit (GPU) computations it is
float32. EchoTorch is not present in the comparison as the version available
on Github was not running on the examples provided. The versions tested are as
follows: PyRCN v0.0.16, ReservoirPy v0.3.2, ReservoirComputing.jl v0.8 release
candidate and pytorch-esn retrieved from Github on March 2022. All the
simulations were run on a Dell XPS 9510 fitted with an Intel i9-11900H CPU, a
Nvidia GeForce RTX 3050 Ti GPU and 16 GB of RAM.

Figure 2: Speed comparison of RC libraries.

| Code quality | Models: ESN | Training | Output | Miscellanea
---|---|---|---|---|---
|  Language |  Documentation |  API/Tutorials |  Tests |  GPU acceleration |  Leaky neurons |  Gated units (Di Sarli et al., 2020) |  Hybrid ESN (Pathak et al., 2018b) |  Multiple reservoirs |  Deep architecture (Gallicchio et al., 2018) |  Non-ESN models |  Linear regression |  Gaussian regression |  SV regression |  Online training |  Generative |  Predictive |  States modifications |  Data sets |  Optimization
ReservoirComputing.jl | Julia | ✓ | ✓/✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✓ | ✗ | ✗
ReservoirPy | Python | ✓ | ✓/✓ | ✓ | ✗ | ✓ | ✗ | ✗ | ✗ | ✓ | ✓ | ✓ | ✗ | ✗ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓
EchoTorch | Python | ✗ | ✓/✗ | ✓ | ✓ | ✓ | ✓ | ✗ | ✓ | ✓ | ✗ | ✓ | ✗ | ✗ | ✗ | ✗ | ✓ | ✗ | ✓ | ✓
pytorch-esn | Python | ✗ | ✓/✗ | ✗ | ✓ | ✓ | ✗ | ✗ | ✗ | ✓ | ✗ | ✓ | ✗ | ✗ | ✓ | ✗ | ✓ | ✗ | ✗ | ✗
PyRCN | Python | ✓ | ✓/✓ | ✓ | ✗ | ✓ | ✗ | ✗ | ✓ | ✓ | ✓ | ✓ | ✗ | ✗ | ✓ | ✗ | ✓ | ✗ | ✓ | ✓

Table 1: Comparison of RC libraries. The _code quality_ section focuses on
high level, quality of life library features. The subsections: _Documentation_
indicates whether the library includes general examples or how-to guides;
_API/Tutorials_ indicates whether API documentation and step by step tutorials
are provided. In the ESN section we list common models included in the
libraries. Subsection _Multiple reservoirs_ indicates whether native multiple
reservoir matrix constructions are provided. In the training section, SV is
short hand for support vector regression. In the _Miscellanea_ section _Data
sets_ indicates whether ready to use data sets are provided in the library;
_Optimization_ stands for native hyperparameter optimization.

##  5 Conclusion and Outlook

ReservoirComputing.jl is a comprehensive, modular and fast library for RC
models with a strong focus on ESNs. In the future we plan to expand the model
library by including LSMs and extreme learning machines (ELMs), while
maintaining the intuitiveness of the current implementation.

## Acknowledgements

This work was partially supported by the German Federal Ministry of Education
and Research (BMBF, 1IS18026B) by funding the competence center for Big Data
and AI “ScaDS.AI” Dresden/Leipzig. K.M. is funded by the iDiv Flexpool
program. M.D.M. and K.M. thank the European Space Agency for funding the
"DeepExtremes" project (AI4Science ITT).

## References

  * Pascanu et al. [2013] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  On the difficulty of training recurrent neural networks.  In _International conference on machine learning_ , pages 1310–1318. PMLR, 2013.
  * Wiggins et al. [2003] Stephen Wiggins, Stephen Wiggins, and Martin Golubitsky.  _Introduction to applied nonlinear dynamical systems and chaos_ , volume 2.  Springer, 2003.
  * Lukoševičius and Jaeger [2009] Mantas Lukoševičius and Herbert Jaeger.  Reservoir computing approaches to recurrent neural network training.  _Computer Science Review_ , 3(3):127–149, 2009\.
  * Jaeger [2001] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  _Bonn, Germany: German National Research Center for Information Technology GMD Technical Report_ , 148(34):13, 2001.
  * Maass et al. [2002] Wolfgang Maass, Thomas Natschläger, and Henry Markram.  Real-time computing without stable states: A new framework for neural computation based on perturbations.  _Neural computation_ , 14(11):2531–2560, 2002\.
  * Gallicchio et al. [2018] Claudio Gallicchio, Alessio Micheli, and Luca Pedrelli.  Design of deep echo state networks.  _Neural Networks_ , 108:33–47, 2018.
  * Chatzis and Demiris [2011] Sotirios P Chatzis and Yiannis Demiris.  Echo state gaussian process.  _IEEE Transactions on Neural Networks_ , 22(9):1435–1445, 2011.
  * Shi and Han [2007] Zhiwei Shi and Min Han.  Support vector echo-state machine for chaotic time-series prediction.  _IEEE transactions on neural networks_ , 18(2):359–372, 2007.
  * Anantharaman et al. [2020] Ranjan Anantharaman, Yingbo Ma, Shashi Gowda, Chris Laughman, Viral Shah, Alan Edelman, and Chris Rackauckas.  Accelerating simulation of stiff nonlinear systems using continuous-time echo state networks.  _arXiv preprint arXiv:2010.04004_ , 2020.
  * Pathak et al. [2018a] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  _Physical review letters_ , 120(2):024102, 2018a.
  * Nadiga [2021] Balasubramanya T Nadiga.  Reservoir computing as a tool for climate predictability studies.  _Journal of Advances in Modeling Earth Systems_ , 13(4):e2020MS002290, 2021.
  * Trouvain et al. [2020] Nathan Trouvain, Luca Pedrelli, Thanh Trung Dinh, and Xavier Hinaut.  Reservoirpy: an efficient and user-friendly library to design echo state networks.  In _International Conference on Artificial Neural Networks_ , pages 494–505. Springer, 2020.
  * Steiner et al. [2021] Peter Steiner, Azarakhsh Jalalvand, Simon Stone, and Peter Birkholz.  Pyrcn: Exploration and application of esns.  _arXiv preprint arXiv:2103.04807_ , 2021.
  * Pontes-Filho et al. [2020] Sidney Pontes-Filho, Pedro Lind, Anis Yazidi, Jianhua Zhang, Hugo Hammer, Gustavo BM Mello, Ioanna Sandvig, Gunnar Tufte, and Stefano Nichele.  A neuro-inspired general framework for the evolution of stochastic dynamical systems: Cellular automata, random boolean networks and echo state networks towards criticality.  _Cognitive Neurodynamics_ , pages 1–18, 2020.
  * Bezanson et al. [2017] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah.  Julia: A fresh approach to numerical computing.  _SIAM review_ , 59(1):65–98, 2017.
  * Konkoli [2018] Zoran Konkoli.  _Reservoir Computing_ , pages 619–629.  Springer US, New York, NY, 2018.  ISBN 978-1-4939-6883-1.  doi:10.1007/978-1-4939-6883-1_683.  URL https://doi.org/10.1007/978-1-4939-6883-1_683.
  * Jaeger [2007] H. Jaeger.  Echo state network.  _Scholarpedia_ , 2(9):2330, 2007.  doi:10.4249/scholarpedia.2330.  revision #196567.
  * Lukoševičius [2012] Mantas Lukoševičius.  A practical guide to applying echo state networks.  In _Neural networks: Tricks of the trade_ , pages 659–686. Springer, 2012.
  * Chattopadhyay et al. [2020] Ashesh Chattopadhyay, Pedram Hassanzadeh, and Devika Subramanian.  Data-driven predictions of a multiscale lorenz 96 chaotic system using machine-learning methods: reservoir computing, artificial neural network, and long short-term memory network.  _Nonlinear Processes in Geophysics_ , 27(3):373–389, 2020.
  * Pathak et al. [2018b] Jaideep Pathak, Alexander Wikner, Rebeckah Fussell, Sarthak Chandra, Brian R Hunt, Michelle Girvan, and Edward Ott.  Hybrid forecasting of chaotic processes: Using machine learning in conjunction with a knowledge-based model.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 28(4):041101, 2018b.
  * Wang et al. [2020] Xinjie Wang, Yaochu Jin, and Kuangrong Hao.  A gated recurrent unit based echo state network.  In _2020 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–7. IEEE, 2020.
  * Lun et al. [2015] Shu-Xian Lun, Xian-Shuang Yao, Hong-Yun Qi, and Hai-Feng Hu.  A novel model of leaky integrator echo state network for time-series prediction.  _Neurocomputing_ , 159:58–66, 2015.
  * Lu et al. [2017] Zhixin Lu, Jaideep Pathak, Brian Hunt, Michelle Girvan, Roger Brockett, and Edward Ott.  Reservoir observers: Model-free inference of unmeasured variables in chaotic systems.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 27(4):041102, 2017.
  * Rodan and Tino [2010] Ali Rodan and Peter Tino.  Minimum complexity echo state network.  _IEEE transactions on neural networks_ , 22(1):131–144, 2010.
  * Rodan and Tiňo [2012] Ali Rodan and Peter Tiňo.  Simple deterministically constructed cycle reservoirs with regular jumps.  _Neural computation_ , 24(7):1822–1852, 2012\.
  * Yang et al. [2018] Cuili Yang, Junfei Qiao, Honggui Han, and Lei Wang.  Design of polynomial echo state networks for time series prediction.  _Neurocomputing_ , 290:148–160, 2018.
  * Yilmaz [2014] Ozgur Yilmaz.  Reservoir computing using cellular automata.  _arXiv preprint arXiv:1410.0162_ , 2014.
  * Nichele and Molund [2017] Stefano Nichele and Andreas Molund.  Deep reservoir computing using cellular automata.  _arXiv preprint arXiv:1703.02806_ , 2017.
  * Martinuzzi [2022] Francesco Martinuzzi.  Martinuzzifrancesco/cellularautomata.jl: v0.0.2, January 2022.  URL https://doi.org/10.5281/zenodo.5879385.
  * Blaom and Vollmer [2020] Anthony D. Blaom and Sebastian J. Vollmer.  Flexible model composition in machine learning and its implementation in MLJ, 2020.
  * Fairbrother et al. [2021] Jamie Fairbrother, Christopher Nemeth, Maxime Rischard, Johanni Brea, and Thomas Pinder.  Gaussianprocesses. jl: A nonparametric bayes package for the julia language.  _Journal of Statistical Software (to appear)_ , 2021.
  * Kornblith and Pastell [2021] Simon Kornblith and Matti Pastell.  https://github.com/juliaml/libsvm.jl, November 2021.  URL https://github.com/JuliaML/LIBSVM.jl.
  * Chang and Lin [2011] Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2:27:1–27:27, 2011.  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * [34] Sciml website.  URL https://sciml.ai/.
  * [35] Sciml at numfocus.  URL https://numfocus.org/project/sciml.
  * Glass and Mackey [2010] L. Glass and M. Mackey.  Mackey-Glass equation.  _Scholarpedia_ , 5(3):6908, 2010.  doi:10.4249/scholarpedia.6908.  revision #186443.
  * Di Sarli et al. [2020] Daniele Di Sarli, Claudio Gallicchio, and Alessio Micheli.  Gated echo state networks: a preliminary study.  In _2020 International Conference on INnovations in Intelligent SysTems and Applications (INISTA)_ , pages 1–5. IEEE, 2020.