This talk is about learning the parameter-to-solution map for systems of partial differential equations (PDEs) that depend on a potentially large number of parameters covering all PDE types for which a {\em stable variational formulation} (SVF) can be found. A central constituent is the notion of {\em variationally correct residual} loss function, meaning that its value is always uniformly proportional to the squared solution error in the norm determined by the SVF, hence facilitating rigorous {\em a posteriori} accuracy control. It is based on a single variational problem, associated with the family of parameter dependent fiber problems, employing the notion of {\em direct integrals of Hilbert spaces}. Since in its original form, the loss function is given as a dual test norm of the residual, a central objective is to develop equivalent computable expressions. A first critical role is played by {\em hybrid hypothesis classes}, whose elements are piecewise polynomial in (low-dimensional) spatio-temporal variables with parameter-dependent coefficients that can be represented, for example, by neural networks. Second, working with first order SVFs, we distinguish two scenarios: (i) the test space can be chosen as an $L_2$-space (such as for elliptic or parabolic problems) so that residuals can be evaluated directly as elements of $L_2$; (ii) when trial and test spaces for the fiber problems depend on the parameters (as for transport equations), we use {\em ultraweak} formulations. In combination with {\em discontinuous Petrov-Galerkin} concepts, the hybrid format is then instrumental to arrive at variationally correct computable residual loss functions. Our findings are illustrated by numerical experiments representing (i) and (ii), namely elliptic boundary value problems with piecewise constant diffusion coefficients and pure transport equations with parameter-dependent convection fields.