Yeah, you could just generate the probabilities of firing with a Linear followed by a Sigmoid, and then call .bernoulli() on the output. That will sample the activations. Then, when you’ll want to backpropagate the errors, you’ll need to provide a reward for every sampled value. Once you do this, you can just call .backward() on the output, and that should do it. You can reuse a single layer multiple times in a single forward pass, it’s perfectly valid in PyTorch.
I think you might have a small bug in your XOR example - you only call the optimizer once after all these iterations. Not sure if that’s what you wanted.