  We develop the first algorithm able to jointly compute the maximum {\it a
posteriori} estimate of the Cosmic Microwave Background (CMB) temperature and
polarization fields, the gravitational potential by which they are lensed, and
cosmological parameters such as the tensor-to-scalar ratio, $r$. This is an
important step towards sampling from the joint posterior probability function
of these quantities, which, assuming Gaussianity of the CMB fields and lensing
potential, contains all available cosmological information and would yield
theoretically optimal constraints. Attaining such optimal constraints will be
crucial for next-generation CMB surveys like CMB-S4, where limits on $r$ could
be improved by factors of a few over currently used sub-optimal quadratic
estimators. The maximization procedure described here depends on a newly
developed lensing algorithm, which we term \textsc{LenseFlow}, and which lenses
a map by solving a system of ordinary differential equations. This description
has conceptual advantages, such as allowing us to give a simple
non-perturbative proof that the lensing determinant is equal to unity in the
weak-lensing regime. The algorithm itself maintains this property even on
pixelized maps, which is crucial for our purposes and unique to
\textsc{LenseFlow} as compared to other lensing algorithms we have tested. It
also has other useful properties such as that it can be trivially inverted
(i.e. delensing) for the same computational cost as the forward operation, and
can be used to compute lensing adjoint, Jacobian, and Hessian operators. We
test and validate the maximization procedure on flat-sky simulations covering
up to 600\,deg$^2$ with non-uniform noise and masking.
