How to use the cvxpy.sum_entries function in cvxpy

# can probably ignore as handled by constraint
    sub_countrol_importance[total_hh_sub_control_index] = 0

    LOG_OVERFLOW = -725
    log_resid_weights = np.log(np.maximum(sub_resid_weights, np.exp(LOG_OVERFLOW))).flatten('F')
    assert not np.isnan(log_resid_weights).any()

    log_parent_resid_weights = \
        np.log(np.maximum(parent_resid_weights, np.exp(LOG_OVERFLOW))).flatten('F')
    assert not np.isnan(log_parent_resid_weights).any()

    # subzone and parent objective and relaxation penalties
    # note: cvxpy overloads * so * in following is matrix multiplication
    objective = cvx.Maximize(
        cvx.sum_entries(cvx.mul_elemwise(log_resid_weights, cvx.vec(x))) +
        cvx.sum_entries(cvx.mul_elemwise(log_parent_resid_weights, cvx.vec(cvx.sum_entries(x, axis=0)))) -  # nopep8
        cvx.sum_entries(relax_le * sub_countrol_importance) -
        cvx.sum_entries(relax_ge * sub_countrol_importance) -
        cvx.sum_entries(cvx.mul_elemwise(parent_countrol_importance, parent_relax_le)) -
        cvx.sum_entries(cvx.mul_elemwise(parent_countrol_importance, parent_relax_ge))
    )

    constraints = [
        (x * sub_incidence) - relax_le >= 0,
        (x * sub_incidence) - relax_le <= lp_right_hand_side,
        (x * sub_incidence) + relax_ge >= lp_right_hand_side,
        (x * sub_incidence) + relax_ge <= hh_constraint_ge_bound,

        x >= 0.0,
        x <= x_max,

        relax_le >= 0.0,

def Hawkes_log_lik(T, alpha_opt, lambda_opt, lambda_ti, survival, for_cvx=False):
    L = 0
    for i in range(len(lambda_ti)):
        if for_cvx and len(lambda_ti) > 0:
            L += CVX.sum_entries(CVX.log(lambda_opt + alpha_opt * lambda_ti[i]))
        else:
            L += np.sum(np.log(lambda_opt + alpha_opt * lambda_ti[i]))

        L -= lambda_opt * T[i] + alpha_opt * survival[i]

    return L

# Get residuals in new marginals from truncating to int
    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
    ]

    prob = cvx.Problem(objective, constraints)

sample_count, control_count = incidence.shape

    # - Decision variables for optimization
    x = cvx.Variable(1, sample_count)

    # - Create positive continuous constraint relaxation variables
    relax_le = cvx.Variable(control_count)
    relax_ge = cvx.Variable(control_count)

    # FIXME - could ignore as handled by constraint?
    control_importance_weights[total_hh_control_index] = 0

    # - Set objective

    objective = cvx.Maximize(
        cvx.sum_entries(cvx.mul_elemwise(log_resid_weights, cvx.vec(x))) -
        cvx.sum_entries(cvx.mul_elemwise(control_importance_weights, relax_le)) -
        cvx.sum_entries(cvx.mul_elemwise(control_importance_weights, relax_ge))
    )

    total_hh_constraint = lp_right_hand_side[total_hh_control_index]

    # 1.0 unless resid_weights is zero
    max_x = (~(resid_weights == 0.0)).astype(float).reshape((1, -1))

    constraints = [
        # - inequality constraints
        cvx.vec(x * incidence) - relax_le >= 0,
        cvx.vec(x * incidence) - relax_le <= lp_right_hand_side,
        cvx.vec(x * incidence) + relax_ge >= lp_right_hand_side,
        cvx.vec(x * incidence) + relax_ge <= hh_constraint_ge_bound,

constraints = [
            x >= 0,
            x.T * hh_table == A,
        ]
        prob = cvx.Problem(objective, constraints)
        prob.solve(solver=cvx.SCS, verbose=verbose_solver)

        return x.value

    else:
        # With relaxation factors
        z = cvx.Variable(n_controls)

        objective = cvx.Maximize(
            cvx.sum_entries(cvx.entr(x) + cvx.mul_elemwise(cvx.log(w.T), x)) +
            cvx.sum_entries(mu * (cvx.entr(z)))
        )

        constraints = [
            x >= 0,
            z >= 0,
            x.T * hh_table == cvx.mul_elemwise(A, z.T),
        ]
        prob = cvx.Problem(objective, constraints)
        prob.solve(solver=cvx.SCS, verbose=verbose_solver)

        return x.value, z.value

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)

"""
            cvxopt.solvers.lp Solves a pair of primal and dual LPs

            minimize    c'*x
            subject to  G*x + s = h
                        A*x = b
                        s >= 0

            maximize    -h'*z - b'*y
            subject to  G'*z + A'*y + c = 0
                        z >= 0.

            cvxopt.solvers.lp(c, G, h[, A, b[, solver[, primalstart[, dualstart]]]])
            """

            obj = cvx.Minimize(cvx.sum_entries(g_ikl * alpha_var))
            prob = cvx.Problem(obj, cons)
            prob.solve()

            alpha_new = np.array(alpha_var.value).T[0]

            """ now the problem collapse into a really simple qp problem """
            # compute beta_new
            for i in range(sample_num):
                alpha_range = class_info.get_range(i)
                alpha_piece = alpha_new[alpha_range[0]:alpha_range[1]]
                c_list = c[i]
                for k in range(class_num):
                    beta_new[k][i] = np.inner(c_list[k], alpha_piece)

            """
            We need to find lambda which will make W(alpha + lambda*alpha_new) be maximum