'N data sets, N param_sets, Polynomial1D'

def test_p1_nset_npset(self):

    p1 = models.Polynomial1D(4, n_models=2)
    y1 = p1(np.array([self.x1, self.x1]).T, model_set_axis=(-1))
    assert (y1.shape == (2, 20))
    assert_allclose(y1[0, :], y1[1, :], atol=(10 ** (-12)))
