Your procedure works but generates a non-uniform density subject to the constraints.
My procedure does neither – in particular the last piece (1-\sum_{j=1}^{k-1}x_j) is not actually constrained.
The Dirichlet distribution has an efficient sampler because of a bijection between the k-simplex and \mathbb{R}^{k-1}. This is not true for a general polytope (which is what it seems you’d like to have).
Honestly, I would use your approach, but swap the uniform distributions with hyper-parametrized betas:
x_i \sim \mathrm{Beta}(\alpha_i, \beta_i) \times \mathrm{min}(\mathrm{resid}, c_i)
\alpha_i \sim \Gamma(a_i^{\alpha}, b_i^{\alpha})
\beta_i \sim \Gamma(a_i^{\beta}, b_i^{\beta})
[something like gamma(2, 2); or you could use an optimizer to make the distribution look uniform] and live with the fact that the prior isn’t uniform.
The most general case: sampling from an arbitrary (convex, hopefully!) polytope; is not trivial; and even estimating volumes is hard. (see http://www.math.cmu.edu/~af1p/Teaching/MCC17/Papers/BF.pdf).
This could be one of those cases where using a potential would work; you could use something like
tt.switch(tt.le(f_constr(x), val), 0, -(val-f_constr(x))*K)
then you could subsequently prune the trace for values outside of the constraint set (since the potentials are imperfect). The alternative is to implement your own interior point method see https://mathoverflow.net/questions/9854/uniformly-sampling-from-convex-polytopes/162327 for a relevant discussion, and a recent paper:
http://www.jmlr.org/papers/volume19/18-158/18-158.pdf
discussing interior point methods for polytope sampling.
Edit:
And using the StickBreaking transformation won’t work either. A simple constraint on the k-simplex winds up being not so simple on the k-1 embedding:
X = sp.stats.dirichlet.rvs(0.45*np.ones(3), 15000)
max_lim = 0.7 * np.ones(3)
within_max = [all(X[i,:] < max_lim) for i in range(X.shape[0])]
Z = forward_val(X)
i = np.where(within_max)[0]
j = np.where(np.logical_not(within_max))[0]
matplotlib.pyplot.scatter(Z[j,0], Z[j,1])
matplotlib.pyplot.scatter(Z[i,0], Z[i,1])
