How to use cvxopt - 10 common examples

re2 = pic.Problem()
XX = re2.add_variable('XX',(6,6),'symmetric')
PP = pic.tools._cplx_mat_to_real_mat(P)
QQ = pic.tools._cplx_mat_to_real_mat(Q)
I = cvx.spdiag([1]*3)
II = pic.tools._cplx_mat_to_real_mat(I)
re2.add_constraint(PP|XX>2)
re2.add_constraint(QQ|XX>2)
re2.add_constraint(XX>>0)
re2.set_objective('min',II|XX)


re3 = pic.Problem()
X = re3.add_variable('X',(3,3),'symmetric')
Y = re3.add_variable('Y',(3,3))
I = cvx.spdiag([1]*3)
re3.add_constraint((P.real()|X)-(P.imag()|Y)>1)
re3.add_constraint((Q.real()|X)-(Q.imag()|Y)>1)
re3.add_constraint(((X & -Y)// (Y & X))>>0)
re3.set_objective('min',I|X)


c2 = pic.Problem()
Z = c2.add_variable('Z',(3,3),'hermitian')
u = c2.add_variable('u',2)
#c2.set_objective('min','I'|Z+3*u[0]+2*u[1])
c2.set_objective('min',Q*Q.H|Z)
c2.add_constraint(P|Z>1)
#c2.add_constraint(Q|Z>1)
#c2.add_constraint(P*Z+Z*P.H+ u[0]*Q*Q.H + u[1]*P*P.H >> P+P.H+Q+Q.H)
c2.add_constraint(R*Z+Z*R.H>>0)
c2.add_constraint('I'|Z==1)

def test_basic_complex(self):
        import cvxopt
        a = cvxopt.matrix([1,-2,3])
        b = cvxopt.matrix([1.0,-2.0,3.0])
        c = cvxopt.matrix([1.0+2j,1-2j,0+1j])
        d = cvxopt.spmatrix([complex(1.0,0.0), complex(0.0,1.0), complex(2.0,-1.0)],[0,1,3],[0,2,3],(4,4))
        e = cvxopt.spmatrix([complex(1.0,0.0), complex(0.0,1.0), complex(2.0,-1.0)],[2,3,3],[1,2,3],(4,4))
        self.assertAlmostEqualLists(list(cvxopt.div(b,c)),[0.2-0.4j,-0.4-0.8j,-3j])
        self.assertAlmostEqualLists(list(cvxopt.div(b,2.0j)),[-0.5j,1j,-1.5j])
        self.assertAlmostEqualLists(list(cvxopt.div(a,c)),[0.2-0.4j,-0.4-0.8j,-3j])
        self.assertAlmostEqualLists(list(cvxopt.div(c,a)),[(1+2j),(-0.5+1j),0.3333333333333333j])
        self.assertAlmostEqualLists(list(cvxopt.div(c,c)),[1.0,1.0,1.0])
        self.assertAlmostEqualLists(list(cvxopt.div(a,2.0j)),[-0.5j,1j,-1.5j])
        self.assertAlmostEqualLists(list(cvxopt.div(c,1.0j)),[2-1j,-2-1j,1+0j])
        self.assertAlmostEqualLists(list(cvxopt.div(1j,c)),[0.4+0.2j,-0.4+0.2j,1+0j])
        self.assertTrue(len(d)+len(e)==len(cvxopt.sparse([d,e])))
        self.assertTrue(len(d)+len(e)==len(cvxopt.sparse([[d],[e]])))

import picos as pic
        P = pic.Problem()
        X = P.add_variable('X',(N,N),'symmetric')
        P.add_constraint(A*X == est*pic.diag_vect(X).T)
        P.add_constraint(A*pic.diag_vect(X)==est)
        P.add_constraint(X<1)
        #P.add_constraint(X>>0)
        P.add_constraint(X>0)
        #P.set_objective('max', (L|X))
        P.set_objective('max', (L|X))
        #P.add_constraint(X[1,4]==0)
        #P.add_constraint(X[2,5]==0)
        #P.add_constraint(X[2,7]==0)
        P.solve()

        H = cvx.sparse([[MMP[:]] for MMP in MP])
        bX = X.value[:]
        Q = pic.Problem()
        lbd = Q.add_variable('lbd',H.size[1],lower=0)
        delta = Q.add_variable('delta',1)
        Q.add_constraint(pic.norm(H*lbd-bX,1)dmax:
                dmax = delta
                Lmax = L
        if delta > 1e-3:
                break


"""
#series of parallel arcs

def compare_opt(self, a_ub, b_ub, c):
        'compare cvx opt versus our glpk interface'

        # make sure we're using floats not ints
        a_ub = [[float(x) for x in row] for row in a_ub]
        b_ub = [float(x) for x in b_ub]
        c = [float(x) for x in c]

        num_vars = len(a_ub[0])

        # solve it with cvxopt
        options = {'show_progress': False}
        sol = cvxopt.solvers.lp(cvxopt.matrix(c), cvxopt.matrix(a_ub).T, cvxopt.matrix(b_ub), options=options)

        if sol['status'] != 'optimal':
            raise RuntimeError("cvxopt LP failed: {}".format(sol['status']))

        res_cvxopt = [float(n) for n in sol['x']]

        # solve it with the glpk <-> hylaa interface
        lp = LpInstance(num_vars, num_vars)

        lp.set_init_constraints(csr_matrix(np.array(a_ub, dtype=float)), np.array(b_ub, dtype=float))
        lp.set_no_output_constraints()
        lp.update_basis_matrix(np.identity(num_vars))

        res_glpk = np.zeros(num_vars)

        lp.minimize(np.array(c, dtype=float), res_glpk)

#     S = A*D^-1*A' + I
        #
        # where D = 2*D1*D2*(D1+D2)^-1, D1 = d[:n]**-2, D2 = d[n:]**-2.

        d1, d2 = W['di'][:n]**2, W['di'][n:]**2

        # ds is square root of diagonal of D
        ds = math.sqrt(2.0) * div( mul( W['di'][:n], W['di'][n:]),
            sqrt(d1+d2) )
        d3 =  div(d2 - d1, d1 + d2)

        # Asc = A*diag(d)^-1/2
        Asc = A * spdiag(ds**-1)

        # S = I + A * D^-1 * A'
        blas.syrk(Asc, S)
        S[::m+1] += 1.0 
        lapack.potrf(S)

        def g(x, y, z):

            x[:n] = 0.5 * ( x[:n] - mul(d3, x[n:]) +
                mul(d1, z[:n] + mul(d3, z[:n])) - mul(d2, z[n:] -
                mul(d3, z[n:])) )
            x[:n] = div( x[:n], ds)

            # Solve
            #
            #     S * v = 0.5 * A * D^-1 * ( bx[:n] -
            #         (D2-D1)*(D1+D2)^-1 * bx[n:] +
            #         D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
            #         D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] )

def test_l1(self):
        setseed(100)
        m,n = 500,250
        P = normal(m,n)
        q = normal(m,1)
        u1,st1 = l1(P,q)
        u2,st2 = l1blas(P,q)
        self.assertTrue(st1 == 'optimal')
        self.assertTrue(st2 == 'optimal')
        self.assertAlmostEqualLists(list(u1),list(u2),places=3)

def test_case3(self):
        m, n = 500, 100
        setseed(100)
        A = normal(m,n)
        b = normal(m)

        x1 = variable(n)
        lp1 = op(max(abs(A*x1-b)))
        lp1.solve()
        self.assertTrue(lp1.status == 'optimal')

        x2 = variable(n)
        lp2 = op(sum(abs(A*x2-b)))
        lp2.solve()
        self.assertTrue(lp2.status == 'optimal')

        x3 = variable(n)
        lp3 = op(sum(max(0, abs(A*x3-b)-0.75, 2*abs(A*x3-b)-2.25)))
        lp3.solve()
        self.assertTrue(lp3.status == 'optimal')

def f(x, y, z):

            # Solve for x[:n]:
            #
            #    A*x[:n] = bx[:n] + P' * ( ((D1-D2)*(D1+D2)^{-1})*bx[n:]
            #        + (2*D1*D2*(D1+D2)^{-1}) * (bz[:m] - bz[m:]) ).

            blas.copy(( mul( div(d1-d2, d1+d2), x[n:]) +
                mul( 2*D, z[:m]-z[m:] ) ), u)
            blas.gemv(P, u, x, beta = 1.0, trans = 'T')
            lapack.potrs(A, x)

            # x[n:] := (D1+D2)^{-1} * (bx[n:] - D1*bz[:m] - D2*bz[m:]
            #     + (D1-D2)*P*x[:n])

            base.gemv(P, x, u)
            x[n:] =  div( x[n:] - mul(d1, z[:m]) - mul(d2, z[m:]) +
                mul(d1-d2, u), d1+d2 )

            # z[:m] := d1[:m] .* ( P*x[:n] - x[n:] - bz[:m])
            # z[m:] := d2[m:] .* (-P*x[:n] - x[n:] - bz[m:])

            z[:m] = mul(di[:m],  u-x[n:]-z[:m])
            z[m:] = mul(di[m:], -u-x[n:]-z[m:])

def fit(self, X, y):
        n_samples, n_features = X.shape

        # Gram matrix
        K = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                K[i,j] = self.kernel(X[i], X[j])

        P = cvxopt.matrix(np.outer(y,y) * K)
        q = cvxopt.matrix(np.ones(n_samples) * -1)
        A = cvxopt.matrix(y, (1,n_samples))
        b = cvxopt.matrix(0.0)

        if self.C is None:
            G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))
            h = cvxopt.matrix(np.zeros(n_samples))
        else:
            tmp1 = np.diag(np.ones(n_samples) * -1)
            tmp2 = np.identity(n_samples)
            G = cvxopt.matrix(np.vstack((tmp1, tmp2)))
            tmp1 = np.zeros(n_samples)
            tmp2 = np.ones(n_samples) * self.C
            h = cvxopt.matrix(np.hstack((tmp1, tmp2)))

        # solve QP problem
        solution = cvxopt.solvers.qp(P, q, G, h, A, b)

def eval_coeff(coeff, params, rows):
    v = coeff.constant_value()
    if v:
        return o.matrix(v, (rows,1))
    else:
        value = 0
        for k,v in coeff.coeff_dict.iteritems():
            if k == '1':
                value += o.matrix(v, (rows,1))
            elif not ismultiply(k):
                p = params[ k.rstrip('\'') ]
                if istranspose(k): value += v*p.T
                else: value += v*p
            else:
                keys = k.split('*')
                keys.reverse()
                mult = 1
                for k1 in keys:
                    p = params[ k.rstrip('\'') ]
                    if istranspose(k): mult = p.T*mult
                    else: mult = p*mult
                value += v*mult
        return value