Correct lithofacies interpretation sourced from wireline log data is an essential source of prior information for joint seismic inversion for facies and impedances, among other applications. However, this information is difficult to interpret or extract manually due to the multivariate and high dimensionality of wireline logs. Facies inference is also challenging for traditional clustering-based approaches because pervasive compaction trends affect a number of petrophysical measurements simultaneously. Another common pitfall in automated clustering approaches is the inability to account for underlying diagenetic processes that correlate with depth. Here, we address these challenges by introducing a rock-physics machine learning toolkit for joint litho-fluid facies classification. The litho-fluid types are inferred from the borehole data within the objective framework of a maximum-likelihood approach for latent facies variables and rock-physics model parameters, explicitly accounting for compaction and depth effects. The inference boils down to an expectation-maximization (EM) algorithm with strong spatial coupling. Each litho-fluid type is associated with an instance of a particular rock-physics model with a unique set of fitting parameters, constrained to a physically reasonable range. These fitting parameters in turn are inferred using bound-constrained optimization as part of the EM algorithm. Outputs produced by the toolkit can be used directly to specify the necessary prior information for seismic inversion, including per-facies rock-physics models and facies proportions. We present an example application of the tool to real borehole data from the North West Shelf of Australia to illustrate the method and discuss its characteristic features in depth.