  Functional forms of the neutron star Equation of State (EoS) are required to
extract the viable EoS band from neutron star observations. Realistic nuclear
EoS, containing deconfined quarks or hyperons, present nontrivial features in
the speed of sound such as bumps, kinks, and plateaus. Using modified Gaussian
processes to model EoS with nontrivial features, we show in a fully Bayesian
analysis incorporating measurements from X-ray sources, gravitational wave
observations, and perturbative QCD results that these features are compatible
with current constraints. We find nontrivial behavior in the EoS plays a role
in understanding the possible phase structure of neutron stars at densities
around 2 $n_{\rm sat}$.
