How to use the cvxpy.mul_elemwise function in cvxpy

def solve_relaxation(prob, *args, **kwargs):
    """Solve the SDP relaxation.
    """

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))
    rel_constr = [X[-1, -1] == 1]

    for f in prob.fs:
        W = f.homogeneous_form()
        lhs = cvx.sum_entries(cvx.mul_elemwise(W, X))
        if f.relop == '==':
            rel_constr.append(lhs == 0)
        else:
            rel_constr.append(lhs <= 0)

    rel_prob = cvx.Problem(rel_obj, rel_constr)
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

    return X.value, rel_prob.value

col_ind += 1

    row = np.concatenate(rows)
    col = np.concatenate(cols)
    data = np.concatenate(datas)
    P = sp.coo_matrix((data, (row, col)), shape=(N, N)).tolil()

    for var in expr.variables():
        col_ind = id_to_col[var.id]
        for j in range(var.size[1]):
            for i in range(var.size[0]):
                P[:, col_ind] -= q
                col_ind += 1

    M = sp.bmat([[P, q], [None, sp.coo_matrix([[r]])]])
    return cvx.sum_entries(cvx.mul_elemwise(M, X))

def inpaint_func(image, mask):
        """Total variation inpainting"""
        inpainted = np.zeros_like(image)
        for c in range(image.shape[2]):
            image_c = image[:, :, c]
            mask_c = mask[:, :, c]
            if np.min(mask_c) > 0:
                # if mask is all ones, no need to inpaint
                inpainted[:, :, c] = image_c
            else:
                h, w = image_c.shape
                inpainted_c_var = cvxpy.Variable(h, w)
                obj = cvxpy.Minimize(cvxpy.tv(inpainted_c_var))
                constraints = [cvxpy.mul_elemwise(mask_c, inpainted_c_var) == cvxpy.mul_elemwise(mask_c, image_c)]
                prob = cvxpy.Problem(obj, constraints)
                # prob.solve(solver=cvxpy.SCS, max_iters=100, eps=1e-2)  # scs solver
                prob.solve()  # default solver
                inpainted[:, :, c] = inpainted_c_var.value
        return inpainted
    return inpaint_func

# TODO: do this efficiently without SDP lifting

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))

    W1 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '<='])
    W2 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '=='])

    rel_prob = cvx.Problem(
        rel_obj,
        [
            cvx.sum_entries(cvx.mul_elemwise(W1, X)) <= 0,
            cvx.sum_entries(cvx.mul_elemwise(W2, X)) == 0,
            X[-1, -1] == 1
        ]
    )
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

    (w, v) = LA.eig(X.value)
    return np.sqrt(np.max(w))*np.asarray(v[:-1, np.argmax(w)]).flatten(), rel_prob.value

Jmat[i, i] = -2
            for i in range(N):
                Jmat[i + 1, i] = 1
                Jmat[i, i + 1] = 1
            Jmat = Jmat / ds ** 2

            # matrix for differentiation
            Amat = np.zeros((N, N + 1))
            for i in range(N):
                Amat[i, i] = -1
                Amat[i, i + 1] = 1
            Amat = Amat / 2 / ds

            xs_var = cvx.Variable(N + 1)
            x_pprime_var = Jmat * xs_var
            s_dddot_var = cvx.mul_elemwise(np.sqrt(xs) + 1e-5, x_pprime_var)
            residue = cvx.Variable()

            constraints = [xs_var >= 0,
                           xs_var[0] == 0,
                           xs_var[-1] == 0,
                           xs_var <= xs,
                           0 <= residue,
                           s_dddot_var - sddd_bnd <= residue,
                           -sddd_bnd - s_dddot_var <= residue]

            obj = cvx.Minimize(1.0 / N * cvx.norm(xs_var - xs, 1)
                               + 1.0 / N * cvx.norm(Amat.dot(xs) - Amat * xs_var, 1)
                               + 1.0 * residue)
            prob = cvx.Problem(obj, constraints)
            status = prob.solve()

# With relaxation factors
    z = cvx.Variable(n_controls, n_tracts)
    q = cvx.Variable(n_controls)

    I = np.ones((n_tracts, 1))

    solved = False
    importance_weights_relaxed = False
    while not solved:
        objective = cvx.Maximize(
            cvx.sum_entries(
                cvx.entr(x) + cvx.mul_elemwise(cvx.log(np.e * w_relative), x)
            ) +
            cvx.sum_entries(
                cvx.mul_elemwise(
                    mu, cvx.entr(z) + cvx.mul_elemwise(cvx.log(np.e), z)
                )
            ) +
            cvx.sum_entries(
                cvx.mul_elemwise(
                    meta_mu, cvx.entr(q) + cvx.mul_elemwise(cvx.log(np.e), q)
                )
            )
        )

        constraints = [
            x >= 0,
            z >= 0,
            q >= 0,
            x * hh_table == cvx.mul_elemwise(A, z.T),
            cvx.mul_elemwise(A.T, z) * I == cvx.mul_elemwise(B.T, q)
        ]

def solve_spectral(prob, *args, **kwargs):
    """Solve the spectral relaxation with lambda = 1.
    """

    # TODO: do this efficiently without SDP lifting

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))

    W1 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '<='])
    W2 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '=='])

    rel_prob = cvx.Problem(
        rel_obj,
        [
            cvx.sum_entries(cvx.mul_elemwise(W1, X)) <= 0,
            cvx.sum_entries(cvx.mul_elemwise(W2, X)) == 0,
            X[-1, -1] == 1
        ]
    )
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

A_residuals = np.dot(x, hh_table) - np.dot(x_int, hh_table)
    x_residuals = x - x_int

    # Coefficients in objective function
    x_log = np.log(x_residuals)

    # Decision variables for optimization
    y = cvx.Variable(n_tracts, n_samples)

    # Relaxation factors
    U = cvx.Variable(n_tracts, n_controls)
    V = cvx.Variable(n_tracts, n_controls)

    objective = cvx.Maximize(
        cvx.sum_entries(
            cvx.sum_entries(cvx.mul_elemwise(x_log, y), axis=1) -
            (gamma) * cvx.sum_entries(U, axis=1) -
            (gamma) * cvx.sum_entries(V, axis=1)
        )
    )

    constraints = [
        y * hh_table <= A_residuals + U,
        y * hh_table >= A_residuals - V,
        U >= 0,
        V >= 0,
        y >= 0,
        y <= 1.0
    ]

    assert not np.isnan(hh_table).any()
    assert not np.isnan(x).any()