code: def getS(self, add_intercept=False):
          
          S = self.S
          if add_intercept is True:
              # S <- cbind(0, rbind(0, S)) # in R
              zeros = np.zeros_like(S[:1, :])
              S = np.vstack([zeros, S])
  
              zeros = np.zeros_like(S[:, :1])
              S = np.hstack([zeros, S])
          return S
docs: Get the penalty matrix S
  
          Returns
              np.array, of shape (n_bases + add_intercept, n_bases + add_intercept)