Removing those does resolve this issue, although I’m not sure why. The function f() does what I was hoping it would with infinite input there. The inf values in gates should disappear outside of the f() function, as they act to allow calculating a probability for all four splits of the input normal distributions across the entire real number range (the erf() handles them appropriately). If I make the values into very small and large numbers (even finfo(float64).eps and finfo(float64).max) it fixes this issue, so that is good.
Ultimately, I am not sure why this fixed the problem because the f() function is returning identical values, I’ll chalk it up to some Theano thing. Thank you for the pointer.