Python cvxpy.Maximize() Examples

def test_entropy_maximization(self):
        set_seed(243)
        n, m, p = 5, 3, 2

        tmp = np.random.rand(n)
        A_np = np.random.randn(m, n)
        b_np = A_np.dot(tmp)
        F_np = np.random.randn(p, n)
        g_np = F_np.dot(tmp) + np.random.rand(p)

        x = cp.Variable(n)
        A = cp.Parameter((m, n))
        b = cp.Parameter(m)
        F = cp.Parameter((p, n))
        g = cp.Parameter(p)
        obj = cp.Maximize(cp.sum(cp.entr(x)) - .01 * cp.sum_squares(x))
        constraints = [A * x == b,
                       F * x <= g]
        prob = cp.Problem(obj, constraints)
        layer = CvxpyLayer(prob, [A, b, F, g], [x])

        A_tch, b_tch, F_tch, g_tch = map(
            lambda x: torch.from_numpy(x).requires_grad_(True), [
                A_np, b_np, F_np, g_np])
        torch.autograd.gradcheck(
            lambda *x: layer(*x, solver_args={"eps": 1e-12,
                                              "max_iters": 10000}),
            (A_tch,
             b_tch,
             F_tch,
             g_tch),
            eps=1e-4,
            atol=1e-3,
            rtol=1e-3) 

def _run_cvx_optimization(self, next_states, rewards, **solver_options):
        """Tensorflow wrapper around a cvxpy value function optimization.

        Parameters
        ----------
        next_states : ndarray
        rewards : ndarray

        Returns
        -------
        values : ndarray
            The optimal values at the states.
        """
        # Define random variables; convert index from np.int64 to regular
        # python int to avoid strange cvxpy error; see:
        # https://github.com/cvxgrp/cvxpy/issues/380
        values = cvxpy.Variable(rewards.shape)

        value_matrix = self.value_function.tri.parameter_derivative(
            next_states)
        # Make cvxpy work with sparse matrices
        value_matrix = cvxpy.Constant(value_matrix)

        objective = cvxpy.Maximize(cvxpy.sum(values))
        constraints = [values <= rewards + self.gamma * value_matrix * values]
        prob = cvxpy.Problem(objective, constraints)

        # Solve optimization problem
        prob.solve(**solver_options)

        # Some error checking
        if not prob.status == cvxpy.OPTIMAL:
            raise OptimizationError('Optimization problem is {}'
                                    .format(prob.status))

        return np.array(values.value) 

def error_bound(L_A, L_b, A_infinity, b_infinity, data, grad_V, eta_jac, decoupling, a_hat, b_hat, controller):
	xs, ts, u_noms, u_perts, V_dot_rs = data
	m = u_noms.shape[1]

	a_hats = [a_hat(x_train, t_train) for x_train, t_train  in zip(xs, ts)]
	b_hats = [b_hat(x_train, t_train) for x_train, t_train  in zip(xs, ts)]
	grad_Vs = [grad_V(x_train, t_train) for x_train, t_train in zip(xs, ts)]
	eta_jacs = [eta_jac(x_train, t_train) for x_train, t_train in zip(xs, ts)]
	V_dot_r_hats = [dot(decoupling(x_train, t_train ), u_pert) + dot(a_hat.T, u_nom + u_pert) + b_hat for x_train, t_train, u_nom, u_pert, a_hat, b_hat in zip(xs, ts, u_noms, u_perts, a_hats, b_hats)]

	def opt(x, t):
		# print(t)

		a = Variable(m)
		b = Variable(1)
		obj = Maximize(a * controller.u(x, t) + b)
		cons = []

		a_hat_test = a_hat(x, t)
		b_hat_test = b_hat(x, t)
		grad_V_test = grad_V(x, t)
		eta_jac_test = eta_jac(x, t)

		def opt_terms(a_hat_train, b_hat_train, grad_V_train, eta_jac_train, u_nom, u_pert, V_dot_r_hat, V_dot_r, x_train):
			u = u_nom + u_pert
			error_obs = abs(V_dot_r - V_dot_r_hat)
			error_model = abs(dot((a_hat_test - a_hat_train).T, u) + b_hat_test - b_hat_train)
			error_inf = norm(dot(grad_V_train, eta_jac_train) - dot(grad_V_test, eta_jac_test))
			error_inf = error_inf * (A_infinity * norm(u) + b_infinity)
			error_lip = min(norm(dot(grad_V_train, eta_jac_train)), norm(dot(grad_V_test, eta_jac_test)))*norm(x_train - x)
			error_lip = error_lip * (L_A * norm(u) + L_b)
			return concatenate([u, ones(1)]), error_obs + error_model + error_inf + error_lip

		zipped = zip(a_hats, b_hats, grad_Vs, eta_jacs, u_noms, u_perts, V_dot_r_hats, V_dot_rs, xs)
		terms = [opt_terms(*params) for params in zipped]
		linear, affine = zip(*terms)
		linear = array(linear)
		linear = concatenate([linear, -linear])
		affine = array(affine)
		affine = concatenate([affine, affine])
		cons = [linear[:, :-1] * a + linear[:, -1] * b <= affine]

		prob = Problem(obj, cons)
		try:
			prob.solve(solver='GLPK', glpk={'msg_lev': 'GLP_MSG_OFF'})
		except Exception:
			print('SOLVER FAILURE', 'State:', x, 'Control:', controller.u(x, t), 'Time:', t)
		return a.value, b.value
	return opt 

def balance_cvx(hh_table, A, w, mu=None, verbose_solver=False):
    """Maximum Entropy allocaion method for a single unit

    Args:
        hh_table (numpy matrix): Table of households categorical data
        A (numpy matrix): Area marginals (controls)
        w (numpy array): Initial household allocation weights
        mu (numpy array): Importance weights of marginals fit accuracy
        verbose_solver (boolean): Provide detailed solver info

    Returns:
        (numpy matrix, numpy matrix): Household weights, relaxation factors
    """

    n_samples, n_controls = hh_table.shape
    x = cvx.Variable(n_samples)

    if mu is None:
        objective = cvx.Maximize(
            cvx.sum_entries(cvx.entr(x) + cvx.mul_elemwise(cvx.log(w.T), x))
        )

        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 diamond_norm_distance(choi0: np.ndarray, choi1: np.ndarray) -> float:
    """
    Return the diamond norm distance between two completely positive
    trace-preserving (CPTP) superoperators, represented as Choi matrices.

    The calculation uses the simplified semidefinite program of Watrous in [CBN]_

    .. note::

        This calculation becomes very slow for 4 or more qubits.

    .. [CBN] Semidefinite programs for completely bounded norms.
          J. Watrous.
          Theory of Computing 5, 11, pp. 217-238 (2009).
          http://theoryofcomputing.org/articles/v005a011
          http://arxiv.org/abs/0901.4709

    :param choi0: A 4**N by 4**N matrix (where N is the number of qubits)
    :param choi1: A 4**N by 4**N matrix (where N is the number of qubits)

    """
    # Kudos: Based on MatLab code written by Marcus P. da Silva
    # (https://github.com/BBN-Q/matlab-diamond-norm/)
    import cvxpy as cvx
    assert choi0.shape == choi1.shape
    assert choi0.shape[0] == choi1.shape[1]
    dim_squared = choi0.shape[0]
    dim = int(np.sqrt(dim_squared))

    delta_choi = choi0 - choi1
    delta_choi = (delta_choi.conj().T + delta_choi) / 2  # Enforce Hermiticity

    # Density matrix must be Hermitian, positive semidefinite, trace 1
    rho = cvx.Variable([dim, dim], complex=True)
    constraints = [rho == rho.H]
    constraints += [rho >> 0]
    constraints += [cvx.trace(rho) == 1]

    # W must be Hermitian, positive semidefinite
    W = cvx.Variable([dim_squared, dim_squared], complex=True)
    constraints += [W == W.H]
    constraints += [W >> 0]

    constraints += [(W - cvx.kron(np.eye(dim), rho)) << 0]

    J = cvx.Parameter([dim_squared, dim_squared], complex=True)
    objective = cvx.Maximize(cvx.real(cvx.trace(J.H * W)))

    prob = cvx.Problem(objective, constraints)

    J.value = delta_choi
    prob.solve()

    dnorm = prob.value * 2

    return dnorm