import numpy as np
import scipy
import cvxpy as cp


def l1_trend_filter(y, vlambda=50, verbose=True):
    """
    Apply l1-trend-filtering to time series y.

    Parameters
    ----------
    vlambda : int, default = 50
        Regularization parameter.
    verbose : bool, default = True
        Whether to print out information.
    """
    n = y.size

    # Form second difference matrix
    e = np.ones((1, n))
    D = scipy.sparse.spdiags(np.vstack((e, -2*e, e)), range(3), n-2, n)

    # Solve l1 trend filtering problem
    x = cp.Variable(shape=n)
    obj = cp.Minimize(0.5 * cp.sum_squares(y - x)
                      + vlambda * cp.norm(D * x, 1))
    prob = cp.Problem(obj)

    # ECOS and SCS solvers fail to converge before the iteration limit. Use CVXOPT instead.
    prob.solve(solver=cp.CVXOPT, verbose=verbose)
    print(f'Solver status: {prob.status}')

    # Check for error.
    if prob.status != cp.OPTIMAL:
        raise Exception('Solver did not converge!')

    print(f'optimal objective value: {obj.value}')

    return x.value