code: def regress_poly(degree, data, remove_mean=True, axis=-1):
      
      timepoints = data.shape[0]
      # Generate design matrix
      X = np.ones((timepoints, 1))  # quick way to calc degree 0
      for i in range(degree):
          polynomial_func = Legendre.basis(i + 1)
          value_array = np.linspace(-1, 1, timepoints)
          X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis]))
      non_constant_regressors = X[:, :-1] if X.shape[1] > 1 else np.array([])
      betas = np.linalg.pinv(X).dot(data)
      if remove_mean:
          datahat = X.dot(betas)
      else:  # disregard the first layer of X, which is degree 0
          datahat = X[:, 1:].dot(betas[1:, ...])
      regressed_data = data - datahat
      return regressed_data, non_constant_regressors
docs: Returns data with degree polynomial regressed out.
      :param bool remove_mean: whether or not demean data (i.e. degree 0),
      :param int axis: numpy array axes along which regression is performed