The current approach without the Jacobian already includes this. The implicit uniform priors on the constrained parameters are already in the posterior p( theta | y ) as defined by the model.
Maybe an example will help. Consider the model:
x ~ exponential(1);
What is the MAP estimate of x? I suspect that you're thinking that it should be at x=0, where the exponential distribution is maximized. However, if you include the Jacobian transformation, you'll get x = 1.
The constrained distribution p( x ) is simple
p(x) = exp(-x)
The unconstrained distribution is
q( x_unc ) = exp(-exp(x_unc)) * exp(x_unc)
= exp(x_unc - exp(x_unc))
If you maximize p(x) wrt x, you get x = 0. If you maximize q(x_unc) wrt x_unc you get x_unc ~= 0 and x = exp(x_unc) ~= 1.
In Stan right now, the optimizer will (try to) give you a value of x=0. (I say try because it would require the optimizer getting to -infinity but that's another issue.)
I don't know where the idea came from that we're not currently doing a valid MAP estimate, we most definitely are.
Perhaps we need to have a skype call to sort this out?