builds the GAM block-diagonal constraint matrix in quadratic form
        out of constraint matrices specified for each feature.

        behaves like a penalty, but with a very large lambda value, ie 1e6.

        Parameters
        ---------
        coefs : array-like containing the coefficients of a term

        constraint_lam : float,
            penalty to impose on the constraint.

            typically this is a very large number.

        constraint_l2 : float,
            loading to improve the numerical conditioning of the constraint
            matrix.

            typically this is a very small number.

        Returns
        -------
        C : sparse CSC matrix containing the model constraints in quadratic form

def build_constraints(self, coef, constraint_lam, constraint_l2):
        """
        builds the GAM block-diagonal constraint matrix in quadratic form
        out of constraint matrices specified for each feature.

        behaves like a penalty, but with a very large lambda value, ie 1e6.

        Parameters
        ---------
        coefs : array-like containing the coefficients of a term

        constraint_lam : float,
            penalty to impose on the constraint.

            typically this is a very large number.

        constraint_l2 : float,
            loading to improve the numerical conditioning of the constraint
            matrix.

            typically this is a very small number.

        Returns
        -------
        C : sparse CSC matrix containing the model constraints in quadratic form
        """
        if self.isintercept:
            return np.array([[0.]])

        Cs = []
        for constraint in self.constraints:

            if constraint is None:
                constraint = 'none'
            if constraint in CONSTRAINTS:
                constraint = CONSTRAINTS[constraint]

            C = constraint(self.n_coefs, coef) * constraint_lam
            Cs.append(C)

        Cs = np.sum(Cs)

        # improve condition
        if Cs.nnz > 0:
            Cs += sp.sparse.diags(constraint_l2 * np.ones(Cs.shape[0]))

        return Cs