'Fit the smoothing spline to a set of (x,y) pairs.
INPUTS:
y       -- response variable
x       -- if None, uses self.x
weights -- optional array of weights
pen     -- constant in front of Gram matrix
OUTPUTS: None
The smoothing spline is determined by self.coef,
subsequent calls of __call__ will be the smoothing spline.
ALGORITHM:
Formally, this solves a minimization:
fhat = ARGMIN_f SUM_i=1^n (y_i-f(x_i))^2 + pen * int f^(2)^2
int is integral. pen is lambda (from Hastie)
See Chapter 5 of
Hastie, Tibshirani and Friedman (2001). "The Elements of Statistical
Learning." Springer-Verlag. 536 pages.
for more details.
TODO:
Should add arbitrary derivative penalty instead of just
second derivative.'

def fit(self, y, x=None, weights=None, pen=0.0):

    banded = True
    if (x is None):
        x = self._x
        bt = self._basisx.copy()
    else:
        bt = self.basis(x)
    if (pen == 0.0):
        banded = False
    if (x.shape != y.shape):
        raise ValueError("x    and    y    shape    do    not    agree,    by    default    x    are                                                                the    Bspline's    internal    knots")
    if (pen >= self.penmax):
        pen = self.penmax
    if (weights is not None):
        self.weights = weights
    else:
        self.weights = 1.0
    _w = np.sqrt(self.weights)
    bt *= _w
    mask = np.flatnonzero((1 - np.alltrue(np.equal(bt, 0), axis=0)))
    bt = bt[:, mask]
    y = y[mask]
    self.df_total = y.shape[0]
    bty = np.squeeze(np.dot(bt, (_w * y)))
    self.N = y.shape[0]
    if (not banded):
        self.btb = np.dot(bt, bt.T)
        _g = _band2array(self.g, lower=1, symmetric=True)
        (self.coef, _, self.rank) = L.lstsq((self.btb + (pen * _g)), bty)[0:3]
        self.rank = min(self.rank, self.btb.shape[0])
        del _g
    else:
        self.btb = np.zeros(self.g.shape, np.float64)
        (nband, nbasis) = self.g.shape
        for i in range(nbasis):
            for k in range(min(nband, (nbasis - i))):
                self.btb[(k, i)] = (bt[i] * bt[(i + k)]).sum()
        bty.shape = (1, bty.shape[0])
        self.pen = pen
        (self.chol, self.coef) = solveh_banded((self.btb + (pen * self.g)), bty, lower=1)
    self.coef = np.squeeze(self.coef)
    self.resid = ((y * self.weights) - np.dot(self.coef, bt))
    self.pen = pen
    del bty
    del mask
    del bt
