code: def transform(self, X):
          
          N, D = X.shape
  
          pow_arr = np.arange(self.order) + 1
  
          # Polynomial terms
          Phi = X[:, :, np.newaxis] ** pow_arr
  
          # Flatten along last axes
          Phi = Phi.reshape(N, D * self.order)
  
          # Prepend intercept
          if self.include_bias:
              Phi = np.hstack((np.ones((N, 1)), Phi))
  
          # TODO: Using np.hstack is about 4x slower than initializing, say,
          # an N by d*order+1 ndarray of ones and assigning the remaining
          # N by d*order values. May want to revisit this implementation.
  
          return Phi
docs: Return this basis applied to X.
  
          Parameters
          ----------
          X: ndarray
              of shape (N, d) of observations where N is the number of samples,
              and d is the dimensionality of X.
  
          Returns
          -------
          ndarray:
              of shape (N, d*order+1), the extra 1 is from a prepended ones
              column.