def _gram_matrix(self):
    knots = self._evaluation_knots()
    c = np.zeros(len(knots))
    ppoly_lst = []
    no_0_intervals = np.where((np.diff(knots) > 0))[0]
    for i in range(self.n_basis):
        c[i] = 1
        pp = PPoly.from_spline((knots, c, (self.order - 1)))
        pp_coefs = pp.c[:, no_0_intervals]
        ppoly_lst.append(pp_coefs)
        c[i] = 0
    matrix = np.zeros((self.n_basis, self.n_basis))
    for interval in range(len(no_0_intervals)):
        for i in range(self.n_basis):
            poly_i = np.trim_zeros(ppoly_lst[i][:, interval], 'f')
            square = polymul(poly_i, poly_i)
            integral = polyint(square)
            matrix[(i, i)] += np.diff(polyval(integral, (self.knots[interval:(interval + 2)] - self.knots[interval])))[0]
            for j in range((i + 1), min((i + self.order), self.n_basis)):
                poly_j = np.trim_zeros(ppoly_lst[j][:, interval], 'f')
                integral = polyint(polymul(poly_i, poly_j))
                matrix[(i, j)] += np.diff(polyval(integral, (self.knots[interval:(interval + 2)] - self.knots[interval])))[0]
                matrix[(j, i)] = matrix[(i, j)]
    return matrix
def _gram_matrix(self):
    knots = self._evaluation_knots()
    c = np.zeros(len(knots))
    ppoly_lst = []
    no_0_intervals = np.where((np.diff(knots) > 0))[0]
    for i in range(self.n_basis):
        c[i] = 1
        pp = PPoly.from_spline((knots, c, (self.order - 1)))
        pp_coefs = pp.c[:, no_0_intervals]
        ppoly_lst.append(pp_coefs)
        c[i] = 0
    matrix = np.zeros((self.n_basis, self.n_basis))
    for interval in range(len(no_0_intervals)):
        for i in range(self.n_basis):
            poly_i = np.trim_zeros(ppoly_lst[i][:, interval], 'f')
            square = polymul(poly_i, poly_i)
            integral = polyint(square)
            matrix[(i, i)] += np.diff(polyval(integral, (self.knots[interval:(interval + 2)] - self.knots[interval])))[0]
            for j in range((i + 1), min((i + self.order), self.n_basis)):
                poly_j = np.trim_zeros(ppoly_lst[j][:, interval], 'f')
                integral = polyint(polymul(poly_i, poly_j))
                matrix[(i, j)] += np.diff(polyval(integral, (self.knots[interval:(interval + 2)] - self.knots[interval])))[0]
                matrix[(j, i)] = matrix[(i, j)]
    return matrix