<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def nmf(Y, A, S, W=None, prox_A=operators.prox_plus, prox_S=operators.prox_plus, proxs_g=None, steps_g=None, Ls=None, slack=0.9, update_order=None, steps_g_update='steps_f', max_iter=1000, e_rel=1e-3, e_abs=0, traceback=None):
 """Non-negative matrix factorization. This method solves the NMF problem minimize || Y - AS ||_2^2 under an arbitrary number of constraints on A and/or S. Args: Y: target matrix MxN A: initial amplitude matrix MxK, will be updated S: initial source matrix KxN, will be updated W: (optional weight matrix MxN) prox_A: direct projection contraint of A prox_S: direct projection constraint of S proxs_g: list of constraints for A or S for ADMM-type optimization steps_g: specific value of step size for proxs_g (experts only!) Ls: list of linear operators for the constraint functions proxs_g If set, needs to have same format as proxs_g. Matrices can be numpy.array, scipy.sparse, or None (for identity). slack: tolerance for (re)evaluation of Lipschitz constants See Steps_AS() for details. update_order: list of factor indices in update order j=0 -> A, j=1 -> S max_iter: maximum iteration number, irrespective of current residuals e_rel: relative error threshold for primal and dual residuals e_abs: absolute error threshold for primal and dual residuals traceback: utils.Traceback to hold variable histories Returns: converged: convence test for A,S errors: difference between latest and previous iterations for A,S See also: algorithms.bsdmm for update_order and steps_g_update utils.AcceleratedProxF for Nesterov acceleration Reference: Moolekamp & Melchior, 2017 (arXiv:1708.09066) """


    # create stepsize callback, needs max of W
    if W is not None:
        # normalize in pixel and band directions to have similar update speeds
        WA = normalizeMatrix(W, 1)
        WS = normalizeMatrix(W, 0)
    else:
        WA = WS = 1
    steps_f = Steps_AS(WA=WA, WS=WS, slack=slack)

    # gradient step, followed by direct application of prox_S or prox_A
    from functools import partial
    f = partial(prox_likelihood, Y=Y, WA=WA, WS=WS, prox_S=prox_S, prox_A=prox_A)

    X = [A, S]
    # use accelerated block-PGM if there's no proxs_g
    if proxs_g is None or not utils.hasNotNone(proxs_g):
        return algorithms.bpgm(X, f, steps_f, accelerated=True, update_order=update_order, max_iter=max_iter, e_rel=e_rel, traceback=traceback)
    else:
        return algorithms.bsdmm(X, f, steps_f, proxs_g, steps_g=steps_g, Ls=Ls, update_order=update_order, steps_g_update=steps_g_update, max_iter=max_iter, e_rel=e_rel, e_abs=e_abs, traceback=traceback)