  Gaussian processes provide a promising framework by which to extrapolate the
equation of state (EoS) of cold, catalyzed matter beyond $1-2$ times nuclear
saturation density. Here we discuss how to extend Gaussian processes to include
nontrivial features in the speed of sound, such as bumps, kinks, and plateaus,
which are predicted by nuclear models with exotic degrees of freedom. Using a
fully Bayesian analysis incorporating measurements from X-ray sources,
gravitational wave observations, and perturbative QCD results, we show that
these features are compatible with current constraints and report on how the
features affect the EoS posteriors.
