Builds a penalty matrix for P-Splines with continuous features.
    Penalizes violation of monotonicity in the feature function.

    Parameters
    ----------
    n : int
        number of splines
    coef : array-like
        coefficients of the feature function
    increasing : bool, default: True
        whether to enforce monotic increasing, or decreasing functions
    Returns
    -------
    penalty matrix : sparse csc matrix of shape (n,n)

def monotonicity_(n, coef, increasing=True):
    """
    Builds a penalty matrix for P-Splines with continuous features.
    Penalizes violation of monotonicity in the feature function.

    Parameters
    ----------
    n : int
        number of splines
    coef : array-like
        coefficients of the feature function
    increasing : bool, default: True
        whether to enforce monotic increasing, or decreasing functions
    Returns
    -------
    penalty matrix : sparse csc matrix of shape (n,n)
    """
    if n != len(coef.ravel()):
        raise ValueError('dimension mismatch: expected n equals len(coef), '\
                         'but found n = {}, coef.shape = {}.'\
                         .format(n, coef.shape))

    if n==1:
        # no monotonic penalty for constant functions
        return sp.sparse.csc_matrix(0.)

    if increasing:
        # only penalize the case where coef_i-1 > coef_i
        mask = sp.sparse.diags((np.diff(coef.ravel()) < 0).astype(float))
    else:
        # only penalize the case where coef_i-1 < coef_i
        mask = sp.sparse.diags((np.diff(coef.ravel()) > 0).astype(float))

    derivative = 1
    D = sparse_diff(sp.sparse.identity(n).tocsc(), n=derivative) * mask
    return D.dot(D.T).tocsc()