With the upcoming generation of galaxy surveys such as Euclid, LSST (Vera Rubin Observatory), DESI, Nancy Roman, and SKAO, among others, the cosmological community will obtain groundbreaking measurements of the large scale structure of the Universe, with a precision never achieved before. This data will contain information about the expansion rate, the growth of non-linear structures, and the presence of standard, and possibly non-standard, species in the matter-energy budget of the Universe. This surge of data promises to unlock the mysteries of Dark Matter and Dark Energy and address the persistent H0 and S8 discrepancies, especially when cross-correlations of different data samples are taken into account. In order to obtain robust results that will survive the scrutiny of time, many previously ignored systematics will have to be taken into account in our theoretical modelling, such as a consistent massive neutrino sector, higher-order galaxy biases, and a careful understanding of baryonic feedback, just to name a few. I will discuss how in the context of the Euclid collaboration, we have implemented these ingredients into the photometric and spectroscopic likelihoods and the tests we are currently performing in the different science working groups. Towards the end I will also address the challenges we face in computational complexity and how this is being resolved with the increased use of emulators and the emerging field of differentiable programming.