Data-driven models, mostly based on statistical and machine learning methods, intend to overcome the limitations of traditional constitutive modelling by directly learning from data. To ensure consistency, recent advancements have introduced physics- and thermodynamics-informed learning frameworks, incorporating first principles through tailored architectures and informative priors. However, the effectiveness of those approaches is hindered by the use of learning biases, where physical constraints are introduced via soft, penalty terms in loss functions, rather than being strictly enforced through hard constraints. This is the case for the second law of thermodynamics, i.e., the entropy production inequality, whose fulfilment by construction remains an open challenge, unlike the first law. One possible solution is to recast the discovery of constitutive equations as the identification of thermodynamic potentials, from which state functions and evolution equations for state variables can be unequivocally derived within the framework of generalised standard materials. While this approach can effectively capture a variety of inelastic material behaviours, the underlying assumption of associative flow rules for the state variables hinders its extension to the broadest range of materials. To address these issues, we develop an approach capable of discovering constitutive equations from scarce and incomplete observations -- a common constraint in traditional experimental observations -- while also satisfying the first and second law of thermodynamics as hard constraints. The capabilities of the method are demonstrated through several representative cases involving in-silico measurement data of inelastic materials, including both associated and non-associated flow rules.