# An empirical study of dust properties at the earliest epochs

Joris Witstok1,2, Gareth C. Jones3, Roberto Maiolino1,2,4, Renske Smit5, and
Raffaella Schneider6,7,8
1Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK
2Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue,
Cambridge CB3 0HE, UK
3Department of Physics, University of Oxford, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK
4Department of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, UK
5Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
6Dipartimento di Fisica, Sapienza, Università di Roma, Piazzale Aldo Moro 5,
00185 Roma, Italy
7INFN, Sezione di Roma I, Piazzale Aldo Moro 2, 00185 Roma, Italy
8INAF/Osservatorio Astronomico di Roma, Via di Frascati 33, 00040 Monte Porzio
Catone, Italy E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted —. Received —; in original form —)

###### Abstract

We present an empirical analysis of the properties of dust-continuum emission
in a sample of 17 galaxies in the early Universe ($4<z<8$) with well-sampled
far-infrared (FIR) spectral energy distributions (SEDs) compiled from the
literature. We place our results into context by self-consistently comparing
to samples of nearby star-forming galaxies, luminous infrared galaxies
(LIRGs), and quasars. With the exception of two sources, we find no
significant evolution in the dust emissivity index across cosmic time,
measuring a consistent value of $\beta_{\text{IR}}=1.8\pm 0.3$ at $z>4$,
suggesting the effective dust properties do not change dramatically for most
galaxies. Despite having comparable stellar masses, we find the high-redshift
galaxies to be similar to, or even more extreme than, LIRGs in the HERUS
sample in terms of dust temperature ($T_{\text{dust}}>40\,\mathrm{K}$) and IR
luminosity ($L_{\text{IR}}>10^{11}\,\mathrm{L_{\odot}}$). We find the dust
temperature evolves mildly towards high redshift, though the LIRGs and quasars
exhibit elevated temperatures indicating a more efficient and/or additional
heating mechanism. Where available, we compare stellar-mass estimates to our
inferred dust masses, whose degeneracy with dust temperature can only be
mitigated with a well-constrained SED. In merely half of the cases the dust
yield may be explained by supernovae alone, with four sources ($44\%$)
significantly exceeding a highly optimistic yield where
$M_{\text{dust}}\approx 0.01M_{*}$. We discuss possible explanations for this
apparent inconsistency and potential observational biases in the measurements
of the dust properties of high-redshift galaxies, including in the current IR-
bright sample.

###### keywords:

galaxies: high-redshift – dark ages, reionization, first stars – methods:
observational – ISM: dust, extinction

††pubyear: 2022††pagerange: An empirical study of dust properties at the
earliest epochs–LABEL:lastpage

##  1 Introduction

Cosmic dust grains are a prominent agent in the physical processes governing
galaxy formation and evolution on the scale of the interstellar medium (ISM).
Dust catalyses the formation of molecules (2017MolAs...9....1W;
2018MNRAS.474.1545C) and the fragmentation of gas clouds (2005ApJ...626..627O;
2006MNRAS.369.1437S), two mechanisms that are essential to star formation.
Furthermore, while the dust mass of a galaxy is negligible compared to its
stellar or gas mass (1998ApJ...501..643D; 2007ApJ...663..866D;
2011A&A...533A..16W), dust grains absorb a significant fraction of the optical
and ultraviolet (UV) light and in the infrared (IR) thermally re-emit the
absorbed energy (e.g. 1979ARA&A..17...73S; 1989ESASP.290...93D;
1999ApJ...521...64M; 2000ApJ...533..682C; 2001ApJ...548..296W;
2003ARA&A..41..241D). This process has been shown to occur even among the
first galaxies to emerge in the Epoch of Reionisation (EoR; e.g.
2017ApJ...851...40L; 2019ApJ...874...27T; 2020MNRAS.493.4294B;
2021Natur.597..489F).

A general consensus has been established on the various channels via which
dust is formed across cosmic time. Specifically, the main sites of dust
creation are thought to be asymptotic giant branch (AGB) stars, supernova (SN)
events, and grain growth in the ISM (e.g. 2015MNRAS.451L..70M;
2020MNRAS.494.1071G; 2022MNRAS.512..989D). However, it is still unclear what
is the exact composition and abundance of dust, particularly in the early
Universe, when the age of the Universe was comparable to typical dust
formation timescales (2001MNRAS.325..726T; 2007MNRAS.378..973B;
2019A&A...624L..13L; 2020MNRAS.497..956S; 2023arXiv230205468W). In addition,
we do not yet fully understand how grain growth occurs in the ISM. It has been
shown that this requires very small grains, with sizes smaller than
$10\,\mathrm{nm}$, in the cold neutral medium (2009ASPC..414..453D;
2016ApJ...831..147Z), and that once the grains are incorporated in dense
molecular clouds, their growth becomes problematic due to the formation of icy
mantles (2016MNRAS.463L.112F; 2018MNRAS.476.1371C). As a result, there exists
significant observational uncertainty on fundamental galaxy properties such as
the star formation rate (SFR), especially when they are inferred exclusively
from rest-frame UV and optical measurements (e.g. 2014ARA&A..52..415M).
Indeed, recent works suggest a significant fraction of obscured star formation
activity may be missed (2020A&A...643A...4F; 2021Natur.597..489F;
2020A&A...643A...3S; 2022MNRAS.512...58F; 2022MNRAS.517.5930S).

With the advent of the Atacama Large Millimeter/submillimeter Array (ALMA),
the first statistical samples of galaxies in the EoR detected by their dust-
continuum emission are starting to be assembled (e.g. 2022ApJ...928...31S;
2022MNRAS.515.3126I). Nearly all of these sources, however, are only observed
in a single photometric band, making it difficult to retrieve properties of
the dust (2021MNRAS.503.4878S; 2022MNRAS.513.3122S). Typically, the spectral
energy distribution (SED) is modelled as a modified blackbody (“greybody”; see
2022MNRAS.515.1751W), assigning ensemble properties to the dust through the
temperature and emissivity parameters, $T_{\text{dust}}$ and
$\beta_{\text{IR}}$, which are commonly assumed for galaxies with few
photometric detections in the far infrared (FIR).

In this work, instead, we compile a list of high-redshift sources with well-
sampled FIR SEDs. We particularly focus on the dust emissivity power-law index
$\beta_{\text{IR}}$, which is connected to the microscopic properties of
grains and therefore potentially holds the key to uncovering their
evolutionary pathways (e.g. 2013A&A...558A..62J). Furthermore, we investigate
the dust mass and temperature, which become degenerate for poorly sampled
SEDs, but have the potential of revealing the main origin of a galaxy’s dust
content built up over relatively short timescales (e.g. 2015A&A...577A..80M;
2020MNRAS.494.1071G; 2022MNRAS.512..989D). In LABEL:sec:Methods, we briefly
describe this sample and two other samples of galaxies in the local Universe.
LABEL:sec:Results_and_discussion describes the results of our analysis and
discusses them in light of dust evolutionary mechanisms and potential changes
in the properties of dust across cosmic time. LABEL:sec:Summary provides the
conclusions of our findings. We adopt the cosmological parameters
$\Omega_{\text{m}}=0.3$, $\Omega_{\Lambda}=0.7$, and
$H_{0}=70\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$ throughout.

Table 1: Properties of local-Universe galaxy samples considered in this work. Sample | Galaxy type | $z$ | $M_{*}\,(\mathrm{M_{\odot}})$ | $A_{\text{dust}}\,(\mathrm{kpc^{2}})$ | $n_{\text{phot}}$ |  References
---|---|---|---|---|---|---
JINGLE | Star-forming | $0.028_{-0.006}^{+0.011}$ | $1.3_{-1.0}^{+2.9}\cdot 10^{10}$ | $1.3_{-0.8}^{+1.9}\cdot 10^{2}$ | $5$ |  2018MNRAS.481.3497S, 2019MNRAS.486.4166S, 2019MNRAS.489.4389L
HERUS | (U)LIRG | $0.086_{-0.043}^{+0.048}$ | $9.0_{-5.0}^{+9.0}\cdot 10^{10}$ | $2.2_{-1.9}^{+2.5}\cdot 10^{3}$ | $5$ |  2003AJ....126.1607S, 2018MNRAS.475.2097C
PG | QSO | $0.14_{-0.08}^{+0.19}$ | – | – | $6$ |  2015ApJS..219...22P
Table 2: Properties of high-redshift sources considered in this work. Source | Type | $z$ | $M_{*}\,(\mathrm{M_{\odot}})$ | $A_{\text{dust}}\,(\mathrm{kpc^{2}})$ | $n_{\text{phot}}$ | References
---|---|---|---|---|---|---
[separator=pipe, late after line= |  |  |  |  |  |
, head to column names]FIR_SED_results.csv  |  |  |  |  |  |


# Adjust factor with volatility model using MAXFLAT low-pass filter and
construct portfolio in China A share market

Ke Zhang
<EMAIL_ADDRESS>
###### Abstract

In the field of quantitative finance, volatility models, such as ARCH, GARCH,
FIGARCH, SV, EWMA, play the key role in risk and portfolio management.
Meanwhile, factor investing is more and more famous since mid of 20 century.
CAPM, Fama–French three-factor model, Fama French five-factor model, MSCI
Barra factor model are mentioned and developed during this period. In this
paper, we will show why we need adjust group of factors by our MAXFLAT low-
pass volatility model. All of our experiments are under China’s CSI 300 and
CSI 500 universe which represent China’s large cap stocks and mid-small cap
stocks. Our result shows adjust factors by MAXFLAT volatility model have
better performance in both large cap and small cap universe than original
factors or other risk adjust factors in China A share. Also the portfolio
constructed by MAXFLAT risk adjust factors have continuous excess return and
lower Beta compare with benchmark index.

Keywords: Volatility model, Multi-factors model, Portfolio construction, China
A share market

##  1 Introduction

Since 1950s, factor model has playing more and more important role in finance
field. Harry Markowitz proposed the Mean-Variance model[9] in 1952 which is
the pioneer of modern finance theory. Then William Sharpe publishes the
capital asset pricing model (CAPM)[10] that firstly introduce the concept of
factor, which explain the portfolio return by excess return and market
exposure. In 1975, Barra mentioned the Barra factor models[21], and
continuously developed it during last 50 years which help us forecast risk for
equity, fixed income, cash and derivative instruments at portfolio and asset
level. In 1993, Eugene Fama and Kenneth French presents the Famam-Fench three
factor model[11] which also include the bond factors. In 2015 Fama extend the
factor model to five factors in his [14]. Factor model are also widely used
and researched in China market[15][16][17][18].

Risk model also play the key role in financial market, in the early of 1980s,
one kind of famous volatility model was introduced by Engle, which is called
Autoregressive conditionally heterosedastic model[2]. Follow the work of
Engle, in 1986, Bollerslev[3] introduced generalized autoregressive
conditional heteroscedasticity model which allowed not only lagged in the
square of errors but also historical variance. Meanwhile, other volatility
model are also developed since then such as stochastic volatility and EWMA
volatility model[7]. Later, people find in any stationary ARCH or GARCH model
that memory decays exponentially fast. Therefore PWMA, IGARCH, FIGARCH[4] or
other long memory volatility were developed to deal with situation has long
memory effect. For example, intra-day high frequency data. However, recently
we can still find sometime people prefer to use rolling truncated period to
estimate volatility which have better performance than all of volatility
models mentioned above. After carefully research, we propose a new volatility
with MAXFLAT loss-pass filter[1][5] that the frequency response of this filter
is flat in the pass-band (i.e. a bandpass filter) and roll-offs towards zero
in the stop-band.. Also widely used rolling period realized volatility is one
of special case of this volatility model. The advantage of this model is that
volatility coefficients will decay slowly at beginning and fast at the end. In
practice, we can use nonlinear state space representation of volatility
model[8][6] with this MAXFLAT filter to estimate our volatility.

In this paper, base on the current research of factors model in China A share
market, we firstly choose and construct several popular factors. Then we do
some data pre-processing work to the stock universe and these factors. And we
use our risk model to adjust these factors. We compare the performance of
factors with different way to adjust. Finally we construct the portfolio with
our risk adjust factors which outperform the benchmark in both large cap and
mid-small cap universe.

##  2 Data, universe and trade assumption

Firstly, our model is low frequency that all of our data is captured at
monthly base and we trade at each end of month. In order to watch the impact
of our model to both large cap stocks and mid-small cap stocks, we pick two
most popular universe in China A share market, CSI 300 and CSI 500. Our
universe is dynamic and point-in time which avoid the survivor bias. On the
other hand, we pick some famous and widely used factors[19][20] from open
source ( https://uqer.datayes.com/). The factors generated from data until
previous day of trading day to avoid look ahead. These factors can be
Categorized as four group, value, profit, growth and tech which can represent
multi-factors model that most researchers and analysts used in China A share
market. We assume commission cost is 5 bps for buy and 1.5 bps for sell. Our
order will totally be filled with 1 cent slippage cost for each side.

Factor | Factor Category
---|---
EP | Value
EB | Value
CFP | Value
CTOP | Value
CFO2EV | Quality
ROE | Quality
ROA | Quality
GrossIncomeRatio | Quality
ARTRate | Quality
DebtsAssetRatio | Quality
OperatingRevenueGrowRate | Growth
NetProfitGrowRate | Growth
SUE | Growth
FEARNG | Growth
FSALESG | Growth
REVS20 | Tech
VOL20 | Tech
ILLIQUIDITY | Tech
Table 1: Common factors we selected

Tab. 1 shows the base style factors we get from open source
https://uqer.datayes.com/.

##  3 Data pre-processing

Each end of month, we get stocks from our universe. However, we will not use
all of these stocks to trade. We remove all the ST and ST* stocks, remove all
the stocks be suspended or non-trade able, remove all the stocks listed less
than 3 month, remove stocks whose 20 average day volume are less 10 million
yuan.

In order to calculate volatility, we need time series of underline asset. We
take log of return series to make it closed to Gaussian distribution.
Meanwhile we adjust the dividend and split for stock return.

As the priority of this paper is not factor analysis or find best factor. We
just use simple rank to weight original factors. Lots of researchers has
mentioned rank factor is a simple but useful and robust way to weight our
portfolio. We firstly winsorize factor to remove outliers, then we rank each
factor and get total score by average of each rank factor. As China A share
market is hard to short, we want to build long only portfolio. Therefore rank
factors automatically give us non-negative weight.

##  4 MAXFLAT low-pass filter

MAXFLAT filter also known as Butterworth filter [1] which is firstly mentioned
by Butterworth at 1930. The MAXFLAT low-pass filter is a type of signal
processing filter designed to have a frequency response that is as flat as
possible in the pass-band and rolls off towards zero in the stop-band. As we
want our volatility Coefficient decay slowly at beginning and fast at end,
MAXFLAT filter easily satisfied our requirement. Also, widely used rolling
window realized volatility is one of spacial case of MAXFLAT filter.

Here is some good property of MAXFLAT low-pass filter.

  * •

Maximal flatness in the pass-band and zero in the stop-band.

  * •

Monotonically decrease from the specified cut-off frequencies.

  * •

Smooth response at all frequencies.

  * •

Half-power frequency that corresponds to the specified cut-off frequencies.

Then let’s see the detail and formula of MAXFLAT low-pass filter.

$|H_{n}(jw)|=\frac{1}{\sqrt{1+w^{2n}}}$ (1)

Eq. 1 is MAXFLAT low-pass frequency response.

Where w is the critical frequency or frequencies and n is the order of the
filter.

Take square and substitute s = jw in Eq. 1, we get

$H(S)H(-S)=\frac{1}{1+{(s/j)}^{2n}}$ (2)

The poles of $H(S)H(-S)$ are given by:

$1+(s/j)^{2n}=0\xrightarrow{}\\\ s^{2n}=-1j^{2n}$ (3)

As

$-1=e^{j\pi{(2k-1)}},j=e^{j\pi/2}$ (4)

We can derive:

$s^{2n}=e^{j\pi(2k-1+n)}$ (5)

Therefore the poles of $H(S)H(-S)$ line on the unity circle:

$s_{k}=e^{\frac{j\pi(2k-1+n)}{2n}}\\\ ,k=1,2,3,...,2n$ (6)

As we are low-pass filter, we only care $H(S)$ which have poles line on the
left side of unity circle. So the poles of our MAXFLAT low-pass filter are:

$s_{k}=e^{\frac{j\pi(2k-1+n)}{2n}}\\\ ,k=1,2,3,...,n$ (7)

Then the transfer function of our low-pass filter is:

$H(s)=\frac{1}{\prod_{i}(s-s_{i})}\\\ ,i=1,2,3,...,n$ (8)

##  5 MAXFLAT low-pass non-linear state space volatility representation

As we already have numerator b and denominator a, we can also get state space
representation in controller canonical form of our filter. As volatility
related to square of return, we need to consider the non-linear state space
representation. To simplify, we use inverse of system representation here.

$\begin{split}Z(t+1)=AZ(t)+By^{2}(t)\\\ x^{2}(t)=y^{2}(t)/CZ(t)\end{split}$
(9)

Eq. 9 is our non-linear State space representation of volatility[8]. Here y(t)
is demean log return at time t (If we didn’t demean y(t) here, there will be a
constant term in state space representation). C,A,B is parameters of state
space model get from MAXFLAT transfer function. x(t) is white noise.

$\displaystyle\sigma^{2}(t)=\sum_{i}{CA^{i}By^{2}(t-i)}$ (10) $\displaystyle
st.CA^{k}B>0,\sum_{k}{CA^{k}B}<1$

Then we can derive square of volatility at time t though Eq. 10

##  6 Volatility model with different decay

Figure 1: Decay of different methods

The above figures shows how fast each methods decay. PWMA decay fastest at
beginning and slowest at the end. EWMA decay faster at beginning than MAXFLAT
and fastest finally. MAXFLAT decay very slow firstly and faster from the
middle.

Coefficient Decay | Volatility Model
---|---
Exponential | EWMA,ARCH,GARCH,SV
Power law | PWMA,IGARCH,FIGARCH
MAXFLAT | Nonlin State Space Vol
Table 2: Vol Model by decay speed

Tab. 2 show the category of some popular volatility model by decay speed.
Recently most of long memory Volatility model are decay similar with power law
speed. Our volatility model with MAXFLAT low-pass filter is a innovation in
volatility field.

##  7 Model construction and comparison

After we pre processed our data, then we consider to construct our portfolio.
Lots of researchers consider volatility as one of factors and linearly add
volatility factors with other factors. We have different opinion that we
linearly weight all Value, Quality, Growth and Tech factors then we divide the
combined factors by our volatility estimator. For volatility model, for
comparison, we pick some famous candidates. Firstly, we choose EWMA who also
known as risk-metric used firstly by JP Morgan. EWMA is also spacial case of
GARCH(1,1). Volatility of EWMA or GARCH model decay exponentially which is too
fast so people consider some long memory volatility model such as PWMA,
FIGARCH[4]. These model mostly decay with power law. In practice, we can also
find lots of researchers or analyst using fix rolling window to estimator
realized volatility. Here we use 1 year and 2 year realized volatility as our
candidates. Finally, we construct a portfolio with adjust factors by MAXFLAT
low-pass filter.

Volatility Methods | Total Return (%) | Sharpe Ratio
---|---|---
EWMA | 120.69 | 0.48
250VOL | 183.48 | 0.61
500VOL | 219.04 | 0.68
PWMA | 243.09 | 0.72
NoVOL | 149.46 | 0.44
MAXFLAT | 242.86 | 0.72
Table 3: Different Volatility Estimator with CSI 300 universe.

In Tab. 3, we construct portfolio with adjust factor by different Volatility
model in universe CSI 300. For each portfolio, we pick first 50 stocks by rank
of adjust factor score. For NoVOL(no Volatility adjustment), it has the worst
Sharpe ratio and second worst total return which tell us we definitely need do
risk adjustment. On the other hand, EWMA result better than NoVol but worst
than other model. Which means, focus on short term volatility is a bad idea
for factor adjustment. MAXFLAT and PWMA has best total return and Sharpe ratio
among all group in this universe.

Volatility Methods | Total Return (%) | Sharpe Ratio
---|---|---
EWMA | 144.13 | 0.45
250VOL | 202.17 | 0.52
500VOL | 206.71 | 0.53
PWMA | 185.49 | 0.50
NoVOL | 151.46 | 0.41
MAXFLAT | 208.51 | 0.53
Table 4: Different Volatility Estimator with CSI 500 universe.

In Tab. 4, here we construct portfolio with adjust factor by different
volatility model in universe CSI 500. For each portfolio, we pick first 50
stocks by rank of adjust factor score. NoVOL(no Volatility adjustment) and
EWMA are still the worst group. Which means, No Volatility adjustment or just
focus on short term volatility is a bad idea for factor adjustment in CSI 500.
MAXFLAT has best total return and Sharpe ratio in all group in CSI 500
universe.

No matter in CSI 300 or CSI 500 universe, MAXFLAT Volatility adjustment show
the best result in term of return and Sharpe ratio. We guess the reason are
volatility clustering effect and long term instead of short term low
volatility of stocks really improve return and Sharpe ratio of portfolio.

##  8 Split MAXFLAT factor portfolio by quantile

In this section, we first construct factor portfolio by rank of MAXFLAT adjust
factor score. Then we split the portfolio to several group by quantile. We
want to see if higher adjust factor score monotonically improve performance of
portfolio.

Figure 2: P$\&$L Quantile with CSI 300 universe using MAXFLAT low-pass
volatility model

In Fig. 2, we construct MAXFLAT factor portfolio in CSI 300. We split the
portfolio to several group by quantile of adjust factor score. In this
universe, we can find group with higher factor score has better total return
monotonically.

Figure 3: P$\&$L Quantile with CSI 500 universe using MAXFLAT volatility model

In Fig. 3, we construct MAXFLAT factor portfolio in CSI 500 universe. In this
universe, except 50 to 100 is slightly worse than 100 to 200, all other group
with higher factor score monotonically has better total return.

As the result above, we showed higher factor score monotonically has better
total return in both large cap and mid-small cap universe in China A share
market.

###  8.1 Compare MAXFLAT factor portfolio with benchmark index

In this section, for CSI 300 and CSI 500 universe, we construct MAXFLAT factor
portfolio with first 50 stocks that have highest MAXFLAT adjust factor score.
Then we compare portfolio performance with corresponding index.

Figure 4: P$\&$L compare with benchmark with CSI 300 universe using MAXFLAT volatility model Porfolio | Total Return (%) | Sharpe | Alpha | Beta
---|---|---|---|---
MAXFLAT | 242.86 | 0.72 | 8.84 | 0.75
CSI 300 Index | 76.94 | 0.29 | 0 | 1
Table 5: Different Volatility Estimator wtih CSI 300 universe.

In Tab. 5, We construct MAXFLAT porfolio and compare with benchmark in CSI 300
universe. MAXFLAT portfolio beat index in term of total return, Sharpe ratio,
Alpha and Beta.

Figure 5: P$\&$L compare with benchmark with CSI 500 universe using MAXFLAT low-pass volatility model Porfolio | Total Return (%) | Sharpe | Alpha | Beta
---|---|---|---|---
MAXFLAT | 208.51 | 0.53 | 7.66 | 0.83
CSI 500 Index | 68.70 | 0.23 | 0 | 1
Table 6: Different Volatility Estimator wtih CSI 500 universe.

In Tab. 6, Here we compare MAXFLAT portfolio with benchmark in CSI 500
universe. MAXFLAT portfolio beat this index for total return, Sharpe ratio,
Alpha and Beta.

##  9 Conclusion

In this paper, we definite a new volatility model by MAXFLAT low-pass filter.
We construct the adjust multi-factors model with different volatility model
for comparison. And we show our MAXFLAT low-pass volatility model outperform
other volatility adjustment in China A share market CSI 300 and CSI 500. The
higher adjust factor score the quantile of portfolio have better performance.
Meanwhile, no matter in universe CSI 300 or CSI 500, our model beat
corresponding benchmark index. The result explain the volatility clustering
effect and the importance of stock long term volatility for portfolio
selection.

## References

  * [1] Butterworth, S. (1930). On the theory of filter amplifiers. Wireless Engineer, 7(6), 536-541.
  * [2] Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the econometric society, 987-1007.
  * [3] Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of econometrics, 31(3), 307-327.
  * [4] Baillie, R. T., Bollerslev, T., Mikkelsen, H. O. (1996). Fractionally integrated generalized autoregressive conditional heteroskedasticity. Journal of econometrics, 74(1), 3-30.
  * [5] Mullhaupt, Andrew P. Triangular Input Balance, Meixner Functions, and a Nearly Rectangular Impulse Response. Unpublished note.
  * [6] Mullhaupt, A. P., Riedel, K. S. (1997). Fast adaptive identification of stable innovation filters. IEEE transactions on signal processing, 45(10), 2616-2619.
  * [7] Higgins, M. L., Bera, A. K. (1992). A class of nonlinear ARCH models. International Economic Review, 137-158.
  * [8] Zhang, Ke (2016). A New Nonlinear State Space Realization of GARCH and Low Rank Adversarial Bandits Algorithm. State University of New York at Stony Brook.
  * [9] Markowitz, H. (1952). Portfolio Selection. The Journal of Finance, 7(1), 77–91.
  * [10] Sharpe, W. F. (1964). Capital asset prices: A theory of market equilibrium under conditions of risk. The journal of finance, 19(3), 425-442.
  * [11] Fama, E. F., French, K. R. (1993). Common risk factors in the returns on stocks and bonds. Journal of financial economics, 33(1), 3-56.
  * [12] Fama, E. F., French, K. R. (1992). The cross‐section of expected stock returns. the Journal of Finance, 47(2), 427-465.
  * [13] Fabozzi, F. J., Markowitz, H. M., Gupta, F. (2008). Portfolio selection. Handbook of finance, 2.
  * [14] Fama, E. F., French, K. R. (2015). A five-factor asset pricing model. Journal of financial economics, 116(1), 1-22.
  * [15] Huang, T. L. (2019). Is the Fama and French five-factor model robust in the Chinese stock market?. Asia Pacific Management Review, 24(3), 278-289.
  * [16] Cai, B. (2022, April). Research on MSCI Barra CNE5 Model and Stock Selection. In 2022 7th International Conference on Social Sciences and Economic Development (ICSSED 2022) (pp. 1505-1512). Atlantis Press.
  * [17] Fang, D., Olteanu-Veerman, D. (2020). The case for factor investing in China a shares. The Journal of Beta Investment Strategies, 11(2), 76-91.
  * [18] Wu, Z., Wang, J., Wang, B. Establishment and Application of Multi-Factor Pricing Model in China A-Shares Market.
  * [19] Zhu, Jiantao (2017). Zhong Mei Shi Chang Yin Zi Xuan Gu Xiao Guo Dui Bi Fen Xi [A comparative analysis of stock factors selection in China and American stock market]. Orient securities.
  * [20] Zhu, Jiantao (2017). Liang Hua Yin Zi Xuan Gu Hui Gu Yu Zhan Wang [Quantitative Factor Stock Selection Review and Prospect]. Orient securities.
  * [21] Sheikh, A. (1996). BARRA’s risk models. Barra Research Insights, 1-24.


# Experiences with managing data parallel computational workflows for High-
throughput Fragment Molecular Orbital (FMO) Calculations

Dimuthu Wannipurage  Indiana UniversityUSA<EMAIL_ADDRESS>,  Indrajit Deb
University of MichiganUSA<EMAIL_ADDRESS>,  Eroma Abeysinghe  Indiana
UniversityUSA<EMAIL_ADDRESS>,  Sudhakar Pamidighantam  Indiana
UniversityBloomingtonIndianaUSA<EMAIL_ADDRESS>,  Suresh Marru  Indiana
UniversityBloomingtonIndianaUSA<EMAIL_ADDRESS>,  Marlon Pierce  Indiana
UniversityUSA<EMAIL_ADDRESS>and  Aaron T. Frank  University of MichiganUSA
<EMAIL_ADDRESS>
(2021)

###### Abstract.

Fragment Molecular Orbital (FMO) calculations provide a framework to speed up
quantum mechanical calculations and so can be used to explore structure-energy
relationships in large and complex biomolecular systems. These calculations
are still onerous, especially when applied to large sets of molecules.
Therefore, cyberinfrastructure that provides mechanisms and user interfaces
that manage job submissions, failed job resubmissions, data retrieval, and
data storage for these calculations are needed. Motivated by the need to
rapidly identify drugs that are likely to bind to targets implicated in SARS-
CoV-2, the virus that causes COVID-19, we developed a static parameter
sweeping framework with Apache Airavata middleware to apply to complexes
formed between SARS-CoV-2 Mpro (the main protease in SARS-CoV-2) and 2820
small-molecules in a drug-repurposing library. Here we describe the
implementation of our framework for managing the executions of the high-
throughput FMO calculations. The approach is general and so should find
utility in large-scale FMO calculations on biomolecular systems.

Covid-19, SARS-CoV-2, Apache Airavata, Science Gateways, FMO Calculations, QM
Calculations

††journalyear: 2021††copyright: acmcopyright††conference: Practice and
Experience in Advanced Research Computing; July 19–22, 2021; Virtual
Conference††booktitle: Practice and Experience in Advanced Research Computing
(PEARC ’21), July 19–22, 2021, Virtual Conference††price: 15.00††doi:
10.1145/3311790.3400848††isbn: 978-1-4503-6689-2/20/07††ccs: General and
reference Experimentation††ccs: Computing methodologies Simulation
environments††ccs: Applied computing Chemistry††ccs: Computing methodologies
Distributed computing methodologies

##  1\. Introduction

SARS-CoV-2 has infected more than 142 million people worldwide and has killed
more than 3 million individuals as of April 2021. Immediately, the repurposing
of old drugs offers us our best chances to reduce the severity of COVID-19 in
patients who have already been infected with SARS-CoV-2.

Figure 1. Our larger vision is to create a scalable Framework for High-
throughput Fragment Molecular Orbital (FMO) Calculations. In this paper we
focus on the data parallel execution aspects.

Scientific goal: We use high-level quantum calculations to guide drug
repurposing efforts for COVID-19. The rationale here is to predict which known
drugs are most likely to bind to known COVID-19 molecular targets. Such drugs
are, in turn, expected to inhibit the replication of SARS-CoV-2 and, in so
doing, stop the spread of COVID-19 (Baker et al., 2020; Bobrowskia et al.,
2020; Filippo, 2020; Duarte et al., 2020; Smith and Smith, 2020). However,
molecular docking methods rely on inaccurate scoring functions that frequently
lead to false-positive and, presumably, false-negative predictions. As an
alternative to the simple empirical scoring functions, quantum mechanical
methods that take into account higher-order quantum mechanical (QM) effects
can be employed. The more accurate QM-derived estimates of drug binding
energies can then be used to re-score and re-prioritize compounds in a given
library based on their docked poses.

Until recently, however, such calculations were prohibitive due to
computational cost. Methodological advances, such as divide-and-conquer
techniques, now make it feasible to calculate drug-binding energies using QM.
In particular, the fragment-molecular-orbital (FMO) method has shown
tremendous promise as a drug screening tool that outperforms approaches
relying on approximate binding energy estimates (Fedorov and Kitaura, 2004;
Okimoto et al., 2018; Heifetz et al., 2020). We, therefore, set out to use the
FMO method to estimate the binding energies of compounds in a drug repurposing
database called SuperDRUG2.0 that contains FDA-approved and marketed drugs
(Siramshetty et al., 2018; Goede et al., 2005; che, 2021). A library of 3993
small-molecule compounds for which 3D structures were available were filtered
to choose small molecules having molecular weight ranging from 200 g/mol to
500 g/mol and atomic numbers more than 10. The library was further filtered to
remove the small-molecules forming salt and in complex with heavy elements.
Finally, a set of 2820 drug-like small molecules was selected and docked to
the active site of SARS-CoV-2 Mpro following standard docking protocol(Ruiz-
Carmona et al., 2014), and the resulting protein-ligand complexes were
subjected to the FMO calculations. This is illustrated as Drug-Repurposing
libary in Fig. 1. Based on the estimates from FMO calculations, the drugs that
are likely to bind to known molecular targets in SARS-CoV-2 can be selected
and advanced as strong candidates for clinical trials.

Computing challenge: The project required the execution of 2820 $\times$ 2 =
5640 independent jobs using the FMO implementation included in the General
Atomic and Molecular Electronic Structure System (GAMESS) (Barca et al., 2020)
an initio quantum chemistry software package. Half of the jobs, corresponding
to the protein-ligand complexes, required 48 processing cores, each with an
execution time of $\sim$7 hrs. The other half, corresponding to the ligand
alone, also used 48 processing cores, each with an execution time of $\sim$5
mins. Job monitoring was required to detect convergence failures, which occur
periodically in FMO calculations. Calculation outputs had to be parsed to
extract energetic and atomic charge data from the GAMESS/FMO log file for
downstream analyses.

By the nature of the execution where the same application is being invoked
with different input configurations, these calculations can be efficiently
performed using a static parameter sweeping technique. In this particular
research use case, input parameters are the various poses of all the above
mentioned FDA approved drug molecules in Cartesian coordinate format, which
are written into text files. These text files become the input parameters for
the GAMESS applications deployed in computational resources to perform FMO
calculations. In this paper we discuss our experiences with developing a
framework (Fig. 1) to submit jobs to computational resources that do not have
an out-of-the-box capability to handle data-parallel parameter sweeping, we
need to have an external job management engine that controls the parameter
sweeping logic. In addition, the framework to manage job executions must be
able to provide a user-friendly gateway interface to submit jobs and monitor
job status, manage job failures and possible re-submissions, and handle large
number of job and post processing final outputs. The framework must also be
generic so that it can be reused for similar and related projects.

##  2\. JOB ORCHESTRATION

We based our implementation on Apache Airavata (Marru et al., 2011), which
provides both a highly flexible job management framework and a turnkey user
interface environment that can be used for managing job executions and data.

Apache Airavata’s job orchestration framework (Wannipurage et al., 2019) is
the component that handles the life cycles of the jobs submitted to High
Performance Computing (HPC) clusters. It consists of multiple workflows
(Figure 2) integrated together to handle the different stages of the life
cycle of a job submitted into a remote computational resource. These job
orchestration workflows are identified as pre-workflows, post-processing
workflows, cancellation workflows and data-parsing workflows. The pre-workflow
is responsible for executing the steps from the job initialization to the job
submission on the computational resource. This includes creating the working
environment on the computational resource for the job, generating the job
script, transferring input files along with the job script and finally
submitting the job to the computational resource. The post-processing workflow
is responsible for staging the output files from the execution host once the
job is completed. A cancellation workflow is executed when user wants to
cancel an already-running job prematurely. A data-parsing workflow runs user
defined data parsers on the output files generated from the job. Typically,
this workflow is executed immediately after the post workflow is completed.
There might be zero or more data-parsing workflows running for a single job
depending on the user or gateway preferences.

Figure 2. Workflow composition of the Job Orchestration Engine of Apache
Airavata

Supporting Static Parameter Sweeping at Job Orchestration: The above job
orchestration framework was designed initially assuming that an experiment
contains a single job. When it comes to static parameter sweeping, there might
be many similar jobs spawned from a single experiment with having changes only
to input parameters. To address this issue, we analyzed how each workflow in
the job orchestrator should be updated to support the required scalability at
the job level of an experiment.

Scaling jobs at the pre-workflow stage: In this particular use case, we have
to submit n jobs in parallel to the target computational cluster with
different input parameters. The pre-workflow submits single parent job to the
cluster, and that job contains the information on how to spin up child jobs
that perform parameter sweeping. The environment setup task contains a batch
of instructions to create the working environments for all those child jobs,
which are executed through a single command in the computational resource.
From the input data side, all the inputs for the parameter sweep are sent as a
tar archive to the remote computational resource, where they are unarchived
and moved to the working environments that were created in the environment
setup phase for each child job.

Application and execution environment: Jobs were run on the Stampede2
supercomputer. We used GAMESS version 2019 R02 and executed using the SKX-
normal partition to utilize the Skylake Intel Xeon Platinum 8160 based nodes
with 48 cores. The GAMESS default executable script, rungms, has been adapted
to perform high throughput parallel processing for multiple job runs, while
each run is assigned a single node. A local scratch directory in Stampede2 is
created that is node-specific for each run. Each run used 24 processor cores
and 24 data servers for processing using MPI parallelization of GAMESS
executable. When a high throughput experiment is submitted in the COVID-QM
gateway in the form of a job array, the nodelist is obtained dynamically from
the SLURM environment variable $SLURM_JOB_NODELIST and used to define the
GAMESS execution host GMS-HOST and a node specific scratch for each of the job
inputs of the job array. Each production job submission contained a job array
of 120 jobs and used the maximum available nodes of 120 per submission
following the reservation allocated for this COVID-19 consortium project.

Our initial approach to scale jobs at the cluster level is to utilize SLURM
job arrays but in Stampede2 where we ran most of the calculations, standard
SLURM job array support was disabled. As a solution, we considered using a
launcher tool that is embedded with the GAMESS application installed in the
cluster with a MPI backend to submit parallel jobs. In this case, the main job
script executed GAMESS child jobs in a loop which runs from 0 to number of
sweep jobs -1 by utilizing the MPI-backed launcher tool. Here, the loop index
is used as the child job index and rest of the configurations were the same as
the proposed job array based solution.

Scaling Post Workflows: Job monitors in the middleware detect when a job
completes, and a post-workflow is executed to pull the generated outputs from
the remote resource into the gateway storage. In a multi-job submission
situation like this, there could be multiple post-workflows running per single
experiment to fetch the outputs of each child job. As the exact job completion
time can not be predicted accurately, post-workflows for these jobs are
triggered randomly. Because of that, there must be a mechanism to guarantee
the eventual consistency of the entire job collection.

In proposed MPI-backed launcher submissions, SLURM only sends email
notifications when the parent job is completed or fails. If the parent job
comes into either of those states, we can safely assume that all the child
jobs also have completed or failed. However, there is no mechanism in SLURM to
send notifications for the status of child jobs. As a solution for this, we
injected a curl command into the beginning and end of the job script with the
child job id as a path parameter. When the job is executed, these curl
commands are also executed. These curl commands point to an HTTP endpoint we
deployed. Using the curl messages received at the HTTP endpoint, we can detect
whether the child job completed or not. This approach enables the opportunity
to fetch the outputs from child jobs before the main job finishes.

Figure 3 depicts the curl message flow and the invocation of the respective
post-workflow for the child job upon receiving the child job completion
notification.

Figure 3. Monitoring of child jobs using injected curl commands to the job
script. Dashed lines and solid lines represent inter-service communications
and data paths, respectively. Numbers on the paths represent the order of
execution.

Handling the eventual experiment completion consistency: Even though curl-
based monitoring helps us find completed child jobs for a parameter sweeping
experiment, it does not give us information about failed jobs. In the case of
a job crash at the computational resource level, we do not get curl messages
as the job script failed before reaching that point. Even in the case of
successful job execution, we can not completely guarantee that curl messages
reaches the HTTP server because of possible firewall restrictions or
intermittent network disruptions. For these reasons, curl-based job monitoring
cannot be used as a mechanism to measure eventual completion consistency.

Email notifications sent by the SLURM scheduler when the parent job is
completed or fails can be used to achieve eventual completion consistency more
reliably but with longer and inconsistent latency. The email monitor of Apache
Airavata fetches this email and passes the status of the parent job into the
job orchestrator to process. As the job orchestrator has information about
already-processed child job post-workflows, it filters out any child job that
was not processed through the curl-based monitoring. Once that subset is
determined, the job orchestrator invokes the post-workflows for those
remaining child jobs. When those are completed, the experiment is marked as
completed or failed based on the parent job status. Figure 4 depicts an
overview of the consistency model described here.

Figure 4.  Eventual consistency model for post workflows.

Output Data Parsing: In some cases, the output of a GAMESS FMO calculation for
certain poses of the drug molecule may not converge. When this occurs, users
should be able to identify these poses and submit jobs with new input
parameters that may allow convergence. We used data parsing workflows to parse
output files and check whether the FMO calculations have converged or not and
also to identify the reason. If the convergence issue is with the coordinates,
we simply energy minimize the systems in our local computer and generate new
input files in which the calculation parameters remain the same, only the
coordinates are different. In case of lack of computational resources, we just
submit a fresh job. The scripts to parse the outputs are packaged as docker
images with defined input and output interfaces. The data-parsing workflow of
Airavata has the capability to load those docker images and run containers
using each job output file as input. After running the designated docker
container, the parsed output for each child job is stored back in the gateway
storage for future analysis.

##  3\. DISCUSSION AND CONCLUSION

The initial testing of the GAMESS application and the pre- and post-parameter
sweep implementations were done using a Jetstream virtual cluster(Coulter et
al., 2019), and construction of the final job scripts as required by the
GAMESS application were configured. The GAMESS output from Stampede2 was
validated against runs on the Jetstream cloud resource and also from a local
computational cluster at the University of Michigan. During testing, it was
determined that using the most recent version of GAMESS was causing
inconsistencies, and GAMESS 2019 R02 was selected as the final version to use
consistently. The GAMESS software was missing atomic radii and parameters for
some of the solute atoms that were required for the calculation of dispersion
and repulsion energies. Those atomic radii and parameters were supplied in the
input explicitly. The computations took about 30 runs with 120 jobs each
including some re-submissions to cover a total 5640 ligand and complex
executions. A small number of jobs ($<$ 0.4%) needed to be rerun as separate
jobs due to convergence issues.

Future directions for this project include:

  1. (1)

Computing the total binding energies and identifying drugs that are predicted
to bind to the SARS-CoV-2 Mpro most favorably.

  2. (2)

Downstream analysis of the pair-wise interaction energies and its
decomposition into energy components to identify key residues at the binding
pocket and the chemical nature of interaction.

  3. (3)

Analyzing the atomic charges to explore ligand-specific polarization effects
at the atomic and molecular level.

  4. (4)

Using the outcome from the above analyses at the lead optimization stage and
also for generating new drug candidates targeting SARS-CoV-2 Mpro.

  5. (5)

Creating a searchable database that stores the key molecular data from GAMESS
log files, e.g., atomic charges, pairwise energies, and solvation free
energies. Access to this data will be important for training machine learning
models capable of reproducing the key results from FMO calculations.

  6. (6)

Using the COVID-QM Gateway to apply FMO calculations to additional
biomolecular systems using proteins, nucleic acids, and their complexes.

###### Acknowledgements.

This work used resources, services, and support provided via the COVID-19 HPC
Consortium (https://covid19-hpc-consortium.org/).

## References

  * (1)
  * che (2021) 2021\.  http://cheminfo.charite.de/superdrug2/index.html.  http://cheminfo.charite.de/superdrug2/index.html
  * Baker et al. (2020) Nancy Baker, Antony Williams, Alexander Tropsha, and Sean Ekins. 2020.  Repurposing Quaternary Ammonium Compounds as Potential Treatments for COVID-19.  (2020).
  * Barca et al. (2020) Giuseppe M. J. Barca, Colleen Bertoni, Laura Carrington, Dipayan Datta, Nuwan De Silva, J. Emiliano Deustua, Dmitri G. Fedorov, Jeffrey R. Gour, Anastasia O. Gunina, Emilie Guidez, Taylor Harville, Stephan Irle, Joe Ivanic, Karol Kowalski, Sarom S. Leang, Hui Li, Wei Li, Jesse J. Lutz, Ilias Magoulas, Joani Mato, Vladimir Mironov, Hiroya Nakata, Buu Q. Pham, Piotr Piecuch, David Poole, Spencer R. Pruitt, Alistair P. Rendell, Luke B. Roskop, Klaus Ruedenberg, Tosaporn Sattasathuchana, Michael W. Schmidt, Jun Shen, Lyudmila Slipchenko, Masha Sosonkina, Vaibhav Sundriyal, Ananta Tiwari, Jorge L. Galvez Vallejo, Bryce Westheimer, Marta Wloch, Peng Xu, Federico Zahariev, and Mark S. Gordon. 2020\.  Recent developments in the general atomic and molecular electronic structure system.  _The Journal of Chemical Physics_ 152, 15 (April 2020), 154102.  https://doi.org/10.1063/5.0005188
  * Bobrowskia et al. (2020) Tesia Bobrowskia, Vinicius M Alvesb, Cleber C Melo-Filhoa, Daniel Korna, Scott Auerbachd, Charles Schmittb, Eugene N Muratova, and Alexander Tropshaa. 2020.  Computational models identify several FDA approved or experimental drugs as putative agents against SARS-CoV-2 2.  (2020).
  * Coulter et al. (2019) J. Eric Coulter, Eroma Abeysinghe, Sudhakar Pamidighantam, and Marlon Pierce. 2019.  Virtual Clusters in the Jetstream Cloud: A Story of Elasticized HPC. In _Proceedings of the Humans in the Loop: Enabling and Facilitating Research on Cloud Computing_ (Chicago, IL, USA) _(HARC ’19)_. Association for Computing Machinery, New York, NY, USA, Article 8, 6 pages.  https://doi.org/10.1145/3355738.3355752
  * Duarte et al. (2020) Rodrigo RR Duarte, Dennis C Copertino Jr, Luis P Iñiguez, Jez L Marston, Douglas F Nixon, and Timothy R Powell. 2020\.  Repurposing FDA-Approved Drugs for COVID-19 Using a Data-Driven Approach.  (2020).
  * Fedorov and Kitaura (2004) Dmitri G Fedorov and Kazuo Kitaura. 2004.  Second order Møller-Plesset perturbation theory based upon the fragment molecular orbital method.  _The Journal of Chemical Physics_ 121, 6 (2004), 2483–2490.
  * Filippo (2020) Juan Di Filippo. 2020\.  In silico Drug Repurposing for COVID-19.  (2020).
  * Goede et al. (2005) Andrean Goede, Mathias Dunkel, Nina Mester, Cornelius Frommel, and Robert Preissner. 2005.  SuperDrug: a conformational drug database.  _Bioinformatics_ 21, 9 (2005), 1751–1753.
  * Heifetz et al. (2020) Alexander Heifetz, Tim James, Michelle Southey, Mike J Bodkin, and Steven Bromidge. 2020.  Guiding Medicinal Chemistry with Fragment Molecular Orbital (FMO) Method.  In _Quantum Mechanics in Drug Discovery_. Springer, 37–48.
  * Marru et al. (2011) Suresh Marru, Lahiru Gunathilake, Chathura Herath, Patanachai Tangchaisin, Marlon Pierce, Chris Mattmann, Raminder Singh, Thilina Gunarathne, Eran Chinthaka, Ross Gardler, Aleksander Slominski, Ate Douma, Srinath Perera, and Sanjiva Weerawarana. 2011.  Apache Airavata: A Framework for Distributed Applications and Computational Workflows. In _Proceedings of the 2011 ACM Workshop on Gateway Computing Environments_ (Seattle, Washington, USA) _(GCE ’11)_. Association for Computing Machinery, New York, NY, USA, 21–28.
  * Okimoto et al. (2018) Noriaki Okimoto, Takao Otsuka, Yoshinori Hirano, and Makoto Taiji. 2018.  Use of the Multilayer Fragment Molecular Orbital Method to Predict the Rank Order of Protein–Ligand Binding Affinities: A Case Study Using Tankyrase 2 Inhibitors.  _ACS Omega_ 3, 4 (2018), 4475–4485.
  * Ruiz-Carmona et al. (2014) Sergio Ruiz-Carmona, Daniel Alvarez-Garcia, Nicolas Foloppe, A Beatriz Garmendia-Doval, Szilveszter Juhos, Peter Schmidtke, Xavier Barril, Roderick E Hubbard, and S David Morley. 2014.  rDock: a fast, versatile and open source program for docking ligands to proteins and nucleic acids.  _PLoS computational biology_ 10, 4 (2014), e1003571.
  * Siramshetty et al. (2018) Vishal B Siramshetty, Oliver Andreas Eckert, Björn-Oliver Gohlke, Andrean Goede, Qiaofeng Chen, Prashanth Devarakonda, Saskia Preissner, and Robert Preissner. 2018.  SuperDRUG2: a one stop resource for approved/marketed drugs.  _Nucleic acids research_ 46, D1 (2018), D1137–D1143.
  * Smith and Smith (2020) Micholas Smith and Jeremy C Smith. 2020.  Repurposing therapeutics for covid-19: Supercomputer-based docking to the sars-cov-2 viral spike protein and viral spike protein-human ace2 interface.  (2020).
  * Wannipurage et al. (2019) Dimuthu Wannipurage, Suresh Marru, Marlon Piece, Eroma Abeysinghe, Sudhakar Pamidighantam, Marcus Christie, Gourav Shenoy, Ajinkya Dhamnaskar, and Lahiru Jayathilaka. 2019.  Implementing a Flexible, Fault Tolerant Job Management System for Science Gateways. In _Proceedings of the Practice and Experience in Advanced Research Computing on Rise of the Machines (Learning)_ (Chicago, IL, USA) _(PEARC ’19)_. Association for Computing Machinery, New York, NY, USA, Article 15, 8 pages.