'Compute Gram inner product matrix, storing it in lower
triangular banded form.
The (i,j) entry is
G_ij = integral b_i^(d) b_j^(d)
where b_i are the basis elements of the BSpline and (d) is the
d-th derivative.
If d is a matrix then, it is assumed to specify a differential
operator as follows: the first row represents the order of derivative
with the second row the coefficient corresponding to that order.
For instance:
[[2, 3],
[3, 1]]
represents 3 * f^(2) + 1 * f^(3).
INPUTS:
d    -- which derivative to apply to each basis element,
if d is a matrix, it is assumed to specify
a differential operator as above
OUTPUTS: gram
gram -- the matrix of inner products of (derivatives)
of the BSpline elements'

def gram(self, d=0):

    d = np.squeeze(d)
    if (np.asarray(d).shape == ()):
        self.g = _hbspline.gram(self.tau, self.m, int(d), int(d))
    else:
        d = np.asarray(d)
        if (d.shape[0] != 2):
            raise ValueError('if    d    is    not    an    integer,    expecting    a    jx2                                                                                array    with    first    row    indicating    order                                                                                of    derivative,    second    row    coefficient    in    front.')
        if (d.shape == (2,)):
            d.shape = (2, 1)
        self.g = 0
        for i in range(d.shape[1]):
            for j in range(d.shape[1]):
                self.g += ((d[(1, i)] * d[(1, j)]) * _hbspline.gram(self.tau, self.m, int(d[(0, i)]), int(d[(0, j)])))
    self.g = self.g.T
    self.d = d
    return np.nan_to_num(self.g)
