Sampling i.i.d. variables with restrictions General Problem: Suppose $X_1,\ldots,X_n \sim \mathbb{P}_X^{\otimes n}$ is a finite sequence of i.i.d. (real- or integer-valued) random variables. Suppose $A\subseteq \mathbb{R}^n$ is a set of "admissible configurations".
Are there efficient methods of sampling from the restriction of $\mathbb{P}_X^{\otimes n}$ onto $A$?
More specific context: I understand that above formulation of the problem is hopelessly general, therefore I would like to give more context on my actual setting. In my case, I am considering a model with $\mathbb{P}_X = \mathrm{Poi}(\lambda)$ and where $A$ is given as the solution set of a finite numbers of equations (representing parity constraints on the Poisson variables). More specifically, I am considering the random current model (cf. this survey article by Duminil-Copin).
Therefore a (perhaps more approachable) version of the problem would be as follows: Consider also a finite set of linear functions $F_1,\ldots,F_N \colon \mathbb{R}^n \to \mathbb{R}$ and an "admissible value set" $V_\mathrm{adm}$. Then let $A:= \{\mathbf{x} \in \mathbb{R}^n\mid F_i(\mathbf{x}) \in V_\mathrm{adm}, \, i=1,\ldots,N\}$.
Is it in this context possible to efficiently sample from the restriction of $\mathbb{P}_X^{\otimes n}$ onto $A$?
I would be glad about input either on the general of more specific formulation of the problem. Also hints to general literature on such problems would be very helpful.
EDIT: As requested, I will try to make the random current setting more explicit: One considers a simple graph $G = (\mathcal{V},\mathcal{E})$. For a subset $S\subseteq \mathcal{V}$, we define the "admissible set" $A_S :=\{ \mathbf{n} = \{\mathbf{n}_e\}_{e\in \mathcal{E}}\in \mathbb{N}_0^\mathcal{E}\mid \mathrm{deg}_\mathbf{n}(v) \text{ odd for $v$ in $S$ and even otherwise}\}$, where $\mathrm{deg}_\mathbf{n}(v) = \sum_{vw \in \mathcal{E}}\mathbf{n}_{vw}$.
 A: A natural approach would be Markov chain Monte Carlo (MCMC), which involves a Markov chain which has your desired measure as its equilibrium measure. At the very least you need this chain to be irreducible, and then you also want it to converge quickly to equilibrium.
If you are content with a sample that is "close to equilibrium" then you could run the chain for a long time, and take the final state as a sample. 
If instead you want do design a way of sampling precisely from equilibrium, it might be a much harder problem, but there are methods such as "coupling from the past" that could have the potential to work. 
Here is an example of the sort of chain that might be suitable for your problem. The analysis of the convergence properties will depend on the structure of your graph. To initialise the chain, you need to find some starting configuration which satisfies all the constraints. Then at each step of the chain, choose, from some suitable distribution, a cycle $v_0 \to v_1 \to \dots \to v_{n-1} \to v_0$ in $G$. Now resample the configuration on the edges of that cycle, according to the conditional distribution of those edges with the rest of the configuration held fixed. (In your model, to keep the parity constraints satisfied, you either need to change the parity of each of the values around the cycle, or to keep the parity of each of the values around the cycle unchanged.) If your cycle is small, then this could just be done by rejection sampling (keep sampling independent Poisson variables until they satsify all the constraints) - or there might also be a more tailored method.
This is some sort of "Glauber dynamics". 
Alternatively, you can consider a more restricted set of transitions (for example, only steps that change the weight of every edge along the cycle by $\pm 1$ - for your model, such changes preserve the required parity constraints), calculating the probability of a move between two given configurations according to the ratio of the equilibrium probabilities of those configurations.  This would be a sort of 
Metropolis-Hastings algorithm. 
EDITED TO ADD:
I realised there is a considerable simplification. You can sample first a model where the state of any edge is "even" or "odd". For this, you again take product measure, conditioned on satisfying the parity constraints at the vertices; now the product measure is the one that puts probability $p$ on "even" and $1-p$ on "odd", where $p$ is the probability that a Poisson($\lambda$) random variable is even. 
Once you have sampled such a configuration (maybe via MCMC) you now assign values to edges, independently. Each "even" edge is sampled from Poisson($\lambda$) conditioned on being even, and each "odd" edge from Poisson($\lambda$) conditioned on being odd. 
For the even/odd model, again a natural move for the MCMC would be to flip all the values around a cycle. More generally, if $x$ is a feasible configuration, then another configuration $y$ is feasible iff their sum (or equivalently difference, since everything is mod 2) $z=x-y$ is feasible for the case where all vertices are even. So a general feasible move is to take such a configuration $z$, and flip all the edges where $z$ is odd. The ratio of probabilities of the states before and after the move is easy to calculate (it's just $p/(1-p)^k$ where $k$ is the net change in the number of even vertices), and that's what you need for Metropolis-Hastings.
(I should stress that none of this actually addresses the question of whether the chain converges to equilibrium in a reasonable time. This may
depend on your graph and on the specifics of how you choose the MCMC updates.)
