I am looking for effective and reliable algorithm which is able to generate random samples of permutations by square doubly stochastic probability matrix $P$ (n x n) distribution ($\sum_{i}p_{i,j} = \sum_{j}p_{i,j} = 1$, $i,j = 1,2,...,n$).
Probability matrix elements $p_{i,j}$ define probability of $i = \pi(j)$ and $0 <= p_{ij} <= 1$.
The sampled permutations should correspond to the defined probability matrix $P$ and probability of each permutation should correspond to the following permutation probability
$$p(\pi)=\frac{\prod_{i} p_{i,\pi(i)}}{\sum_{all \pi}\prod_{j} p_{j,\pi(j)}}=\frac{\prod_{i} p_{i,\pi(i)}}{perm(P)}$$
Example: Probability distribution matrix
$$ P = \left( \begin{array}{ccc} 0 & 1 & 0 \\ 1/2 & 0 & 1/2 \\ 1/2 & 0 & 1/2 \end{array} \right) $$
produce only two permutations $\pi_1 = (2 1 3)$ and $\pi_2 =(2 3 1)$ with the same probability $p(\pi_1) = p(\pi_2) =1/2$
This question is a follow-up to Random permutations by probability matrix
I tried a few variations of node histogram based sampling algorithm (NHBSA) by Shigeyoshi Tsutsui. This method is used as a main tool on permutation optimization problems. But the results are not very promising. So, I ask this question here, to learn what is current state of art on this topic.
I tried to test a few algorithms, but all these methods produce permutations with more or less distorted probability matrix even for small n (3 <= n <= 10 ). This is the reason why I am asking for help.
Some comments here states that interpretation of matrix $P$ and permutations probability $p(\pi)$ are not clear enough to define unique permutations distribution.
So, is possible to use $\prod_{i} p_{i,\pi(i)}$ for ranking of permutation probabilities, where permutations $\pi$ are generated randomly only by (0,1) restriction matrix via Hungarian algorithm for bipartite graph matching?
Comment: I am convinced now, that only "simple" way how to solve my problem is to generate arbitrary (random) permutations $\pi$ which satisfy to (0,1) restrictions in probability matrix only (elements $p_{i,j}>0$ should be transformed to $p_{i,j}=1$) and then rank all these permutations by $\prod_{i} p_{i,\pi(i)}$ probability.
The problem how to effectively generate (0,1) restricted permutations (permutations with restricted positions) is different problem.
