Accurate knowledge of phase coexistence regions, i.e., solubility gaps (SGs), is key to the development of mixed transition metal oxides for various applications, such as thermochemical energy storage, or catalysis. However, predicting a SG from first principles in these materials is particularly challenging due to the complex interplay between several sources of entropy, the large configuration space, and the computational expense of ab initio calculations. We present an approach that yields an accurate prediction of the experimental Hausmannite-spinel SG in the case of (CoxMn1-x)3O4. The method uses machine learning to extend an ab initio dataset of hundreds of structures, and it includes many different entropic contributions to the free energy. We demonstrate and quantify the crucial roles of phonon and paramagnetic entropy, and the importance of sampling higher-energy configurations, and correcting for finite-size limitations in the ab initio supercell configurations.