  Many astrophysical applications require efficient yet reliable forecasts of
stellar evolution tracks. One example is population synthesis, which generates
forward predictions of models for comparison with observations. The majority of
state-of-the-art rapid population synthesis methods are based on analytic
fitting formulae to stellar evolution tracks that are computationally cheap to
sample statistically over a continuous parameter range. The computational costs
of running detailed stellar evolution codes, such as MESA, over wide and
densely sampled parameter grids are prohibitive, while stellar-age based
interpolation in-between sparsely sampled grid points leads to intolerably
large systematic prediction errors. In this work, we provide two solutions for
automated interpolation methods that offer satisfactory trade-off points
between cost-efficiency and accuracy. We construct a timescale-adapted
evolutionary coordinate and use it in a two-step interpolation scheme that
traces the evolution of stars from ZAMS all the way to the end of core helium
burning while covering a mass range from ${0.65}$ to $300 \, \mathrm{M_\odot}$.
The feedforward neural network regression model (first solution) that we train
to predict stellar surface variables can make millions of predictions,
sufficiently accurate over the entire parameter space, within tens of seconds
on a 4-core CPU. The hierarchical nearest-neighbor interpolation algorithm
(second solution) that we hard-code to the same end achieves even higher
predictive accuracy, the same algorithm remains applicable to all stellar
variables evolved over time, but it is two orders of magnitude slower. Our
methodological framework is demonstrated to work on the MIST (Choi et al. 2016)
data set. Finally, we discuss the prospective applications of these methods and
provide guidelines for generalizing them to higher dimensional parameter
spaces.
