def polynomial_features(X, degree):
    (n_samples, n_features) = np.shape(X)

    def index_combinations():
        combs = [combinations_with_replacement(range(n_features), i) for i in range(0, (degree + 1))]
        flat_combs = [item for sublist in combs for item in sublist]
        return flat_combs
    combinations = index_combinations()
    n_output_features = len(combinations)
    X_new = np.empty((n_samples, n_output_features))
    for (i, index_combs) in enumerate(combinations):
        X_new[:, i] = np.prod(X[:, index_combs], axis=1)
    return X_new
def polynomial_features(X, degree):
    (n_samples, n_features) = np.shape(X)

    def index_combinations():
        combs = [combinations_with_replacement(range(n_features), i) for i in range(0, (degree + 1))]
        flat_combs = [item for sublist in combs for item in sublist]
        return flat_combs
    combinations = index_combinations()
    n_output_features = len(combinations)
    X_new = np.empty((n_samples, n_output_features))
    for (i, index_combs) in enumerate(combinations):
        X_new[:, i] = np.prod(X[:, index_combs], axis=1)
    return X_new