  We present a semi-analytic model for predicting kilonova light curves from
the mergers of neutron stars with black holes (NSBH). The model is integrated
into the MOSFiT platform, and can generate light curves from input binary
properties and nuclear equation-of-state considerations, or incorporate
measurements from gravitational wave (GW) detectors to perform multi-messenger
parameter estimation. The rapid framework enables the generation of NSBH
kilonova distributions from binary populations, light curve predictions from GW
data, and statistically meaningful comparisons with an equivalent BNS model in
MOSFiT. We investigate a sample of kilonova candidates associated with
cosmological short gamma-ray bursts, and demonstrate that they are broadly
consistent with being driven by NSBH systems, though most have limited data. We
also perform fits to the very well sampled GW170817, and show that the
inability of an NSBH merger to produce lanthanide-poor ejecta results in a
significant underestimate of the early (< 2 days) optical emission. Our model
indicates that NSBH-driven kilonovae may peak up to a week after merger at
optical wavelengths for some observer angles. This demonstrates the need for
early coverage of emergent kilonovae in cases where the GW signal is either
ambiguous or absent; they likely cannot be distinguished from BNS mergers by
the light curves alone from ~2 days after the merger. We also discuss the
detectability of our model kilonovae with the Vera C. Rubin Observatory's
Legacy Survey of Space and Time (LSST).
