[MRG+1] compute poly features directly

Here's a much simpler alternative to #3194 for fixing #3191. `itertools` and NumPy have all the required functionality.

The great thing about this approach is that we can turn the actual polynomials (x², y²) off and get only the interaction terms (xy) by substituting `combinations` for `combinations_with_replacement`.