Python cvxopt.spmatrix() Examples

def cvxopt_matrix(M):
    """
    Convert matrix to CVXOPT format.

    Parameters
    ----------
    M : numpy.array
        Matrix in NumPy format.

    Returns
    -------
    N : cvxopt.matrix
        Matrix in CVXOPT format.
    """
    if type(M) is ndarray:
        return matrix(M)
    elif type(M) is spmatrix or type(M) is matrix:
        return M
    coo = M.tocoo()
    return spmatrix(
        coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=M.shape) 

def _get_projector(R, N_ex):  # !
    # Projection
    dim = N_ex
    P = spmatrix(1, range(dim), range(dim))
    glast = matrix(np.ones((1, dim)))
    G = sparse([-P, P, glast])
    h1 = np.zeros(dim)
    h2 = np.ones(dim)
    h = matrix(np.concatenate([h1, h2, [R]]))

    def _project(pt):
        print('start projection')
        # pt = gamma.eval()
        q = matrix(- np.array(pt, dtype=np.float64))
        # if np.linalg.norm(pt, ord=1) < R:
        #    return
        _res = cvxopt.solvers.qp(P, q, G, h, initvals=q)
        _resx = np.array(_res['x'], dtype=np.float32)[:, 0]
        # gamma_assign.eval(feed_dict={grad_hyper: _resx})
        return _resx

    return _project


# TODO check the following functions (look right) 

def __init__(self, fun, Exp, funstring):
        Expression.__init__(self, self.funstring + '( ' + Exp.string + ')')
        self.fun = fun
        r"""The function ``f`` applied to the affine expression.
                This function must take in argument a
                :func:`cvxopt sparse matrix <cvxopt:cvxopt.spmatrix>` ``X``.
                Moreover, the call ``fx,grad,hess=f(X)``
                must return the function value :math:`f(X)` in ``fx``,
                the gradient :math:`\nabla f(X)` in the
                :func:`cvxopt matrix <cvxopt:cvxopt.matrix>` ``grad``,
                and the Hessian :math:`\nabla^2 f(X)` in the
                :func:`cvxopt sparse matrix <cvxopt:cvxopt.spmatrix>` ``hess``.
                """
        self.Exp = Exp
        """The affine expression to which the function is applied"""
        self.funstring = funstring
        """a string representation of the function name"""
        #self.string=self.funstring+'( '+self.Exp.affstring()+' )' 

def _break_cols(mat, sizes):
    n = len(sizes)
    I, J, V = [], [], []
    for i in range(n):
        I.append([])
        J.append([])
        V.append([])
    cumsz = np.cumsum(sizes)
    import bisect
    for i, j, v in zip(mat.I, mat.J, mat.V):
        block = bisect.bisect(cumsz, j)
        I[block].append(i)
        V[block].append(v)
        if block == 0:
            J[block].append(j)
        else:
            J[block].append(j - cumsz[block - 1])
    return [spmatrix(V[k], I[k], J[k], (mat.size[0], sz))
            for k, sz in enumerate(sizes)] 

def _break_rows(mat, sizes):
    n = len(sizes)
    I, J, V = [], [], []
    for i in range(n):
        I.append([])
        J.append([])
        V.append([])
    cumsz = np.cumsum(sizes)
    import bisect
    for i, j, v in zip(mat.I, mat.J, mat.V):
        block = bisect.bisect(cumsz, i)
        J[block].append(j)
        V[block].append(v)
        if block == 0:
            I[block].append(i)
        else:
            I[block].append(i - cumsz[block - 1])
    return [spmatrix(V[k], I[k], J[k], (sz, mat.size[1]))
            for k, sz in enumerate(sizes)] 

def _blocdiag(X, n):
    """
    makes diagonal blocs of X, for indices in [sub1,sub2[
    n indicates the total number of blocks (horizontally)
    """
    if not isinstance(X, cvx.base.spmatrix):
        X = cvx.sparse(X)
    if n==1:
        return X
    else:
        Z = spmatrix([],[],[],X.size)
        mat = []
        for i in range(n):
            col = [Z]*(n-1)
            col.insert(i,X)
            mat.append(col)
        return cvx.sparse(mat) 

def diag(self, dim):
        if self.size != (1, 1):
            raise Exception('not implemented')
        selfcopy = self.copy()
        idx = cvx.spdiag([1.] * dim)[:].I

        for k in self.factors:
            tc = 'z' if self.factors[k].typecode=='z' else 'd'
            selfcopy.factors[k] = spmatrix(
                [], [], [], (dim**2, self.factors[k].size[1]),tc=tc)
            for i in idx:
                selfcopy.factors[k][i, :] = self.factors[k]
        if not self.constant is None:
            tc = 'z' if self.constant.typecode=='z' else 'd'
            selfcopy.constant = cvx.matrix(0., (dim**2, 1),tc=tc)
            for i in idx:
                selfcopy.constant[i] = self.constant[0]
        else:
            selfcopy.constant = None
        selfcopy._size = (dim, dim)
        selfcopy.string = 'diag(' + selfcopy.string + ')'
        return selfcopy 

def inplace_transpose(self):
        if isinstance(self, Variable):
            raise Exception(
                'inplace_transpose should not be called on a Variable object')
        for k in self.factors:
            bsize = self.size[0]
            bsize2 = self.size[1]
            I0 = [(i // bsize) + (i % bsize) *
                  bsize2 for i in self.factors[k].I]
            J = self.factors[k].J
            V = self.factors[k].V
            self.factors[k] = spmatrix(V, I0, J, self.factors[k].size)

        if not (self.constant is None):
            self.constant = cvx.matrix(self.constant,
                                       self.size).T[:]
        self._size = (self.size[1], self.size[0])
        if (('*' in self.affstring()) or ('/' in self.affstring())
                or ('+' in self.affstring()) or ('-' in self.affstring())):
            self.string = '( ' + self.string + ' ).T'
        else:
            self.string += '.T' 

def compute_displacements(self, xPhys):
        # Setup and solve FE problem
        sK = ((self.KE.flatten()[np.newaxis]).T * (
            self.Emin + (xPhys)**self.penal *
            (self.Emax - self.Emin))).flatten(order='F')
        K = scipy.sparse.coo_matrix((sK, (self.iK, self.jK)),
            shape=(self.ndof, self.ndof)).tocsc()
        # Remove constrained dofs from matrix and convert to coo
        K = deleterowcol(K, self.fixed, self.fixed).tocoo()
        # Solve system
        K1 = cvxopt.spmatrix(K.data, K.row.astype(np.int), K.col.astype(np.int))
        B = cvxopt.matrix(self.f[self.free, :])
        cvxopt.cholmod.linsolve(K1, B)
        self.u[self.free, :] = np.array(B)[:, :] 

def klu_linsolve(A, b):
    """
    KLU wrapper function for linear system solve A x = b
    :param A: System matrix
    :param b: right hand side
    :return: solution
    """
    A2 = A.tocoo()
    A_cvxopt = cvxopt.spmatrix(A2.data, A2.row, A2.col, A2.shape, 'd')
    x = cvxopt.matrix(b)
    klu.linsolve(A_cvxopt, x)
    return np.array(x)[:, 0] 

def __init__(self, dim, n_labels, mode=1):
        self.mode = mode
        self.dim = dim
        self.n_labels = n_labels
        if mode == 1:
            self.P = spmatrix(1, range(dim * n_labels), range(dim * n_labels))

            glast = matrix(np.ones((1, dim * n_labels)))
            self.G = sparse([-self.P, self.P, glast])

            h1 = np.zeros(dim * n_labels)
            h2 = np.ones(dim * n_labels)
            self.h = matrix(np.concatenate([h1, h2, [dim]]))
        elif mode == 2:
            self.P = spmatrix(1, range(n_labels), range(n_labels))
            glast = matrix(np.ones((1, n_labels)))
            self.G = sparse([-self.P, self.P, glast])

            h1 = np.zeros(n_labels)
            h2 = np.ones(n_labels)
            self.h = matrix(np.concatenate([h1, h2, [1]]))
        elif mode == 3:
            self.P = spmatrix(1, range(n_labels), range(n_labels))
            self.A = matrix(np.ones((1, n_labels)))
            self.G = sparse([-self.P, self.P])

            h1 = np.zeros(n_labels)
            h2 = np.ones(n_labels)
            self.h = matrix(np.concatenate([h1, h2]))
            self.b = matrix(np.ones(1)) 

def inplace_Htranspose(self):
        if isinstance(self, Variable):
            raise Exception(
                'inplace_transpose should not be called on a Variable object')
        for k in self.factors:
            if k.vtype == 'hermitian':
                bsize = self.size[0]
                bsize2 = self.size[1]
                vsize = k.size[0]
                vsize2 = k.size[1]
                I0 = [(i // bsize) + (i % bsize) *
                      bsize2 for i in self.factors[k].I]
                J0 = [(j // vsize) + (j % vsize) *
                      vsize2 for j in self.factors[k].J]
                V0 = [v.conjugate() for v in self.factors[k].V]
                self.factors[k] = spmatrix(
                    V0, I0, J0, self.factors[k].size)
            else:
                F = self.factors[k]
                bsize = self.size[0]
                bsize2 = self.size[1]
                I0 = [(i // bsize) + (i % bsize) * bsize2 for i in F.I]
                J = F.J
                V = [v.conjugate() for v in F.V]
                self.factors[k] = spmatrix(V, I0, J, F.size)

        if not (self.constant is None):
            if self.constant.typecode == 'z':
                self.constant = cvx.matrix(self.constant,
                                           self.size).H[:]
            else:
                self.constant = cvx.matrix(self.constant,
                                           self.size).T[:]
        self._size = (self.size[1], self.size[0])
        if (('*' in self.affstring()) or ('/' in self.affstring())
                or ('+' in self.affstring()) or ('-' in self.affstring())):
            self.string = '( ' + self.string + ' ).H'
        else:
            self.string += '.H' 

def inplace_conjugate(self):
        if isinstance(self, Variable):
            raise Exception(
                'inplace_conjugate should not be called on a Variable object')
        for k in self.factors:
            if k.vtype == 'hermitian':
                fack = self.factors[k]
                I = fack.I
                J = fack.J
                V = fack.V
                J0 = []
                V0 = []
                n = int(fack.size[1]**0.5)
                for j, v in zip(J, V):
                    col, row = divmod(j, n)
                    J0.append(row * n + col)
                    V0.append(v.conjugate())
                self.factors[k] = spmatrix(V0, I, J0, fack.size)

            elif self.factors[k].typecode == 'z':
                Fr = self.factors[k].real()
                Fi = self.factors[k].imag()
                self.factors[k] = Fr - 1j * Fi
        if not self.constant is None:
            if self.constant.typecode == 'z':
                Fr = self.constant.real()
                Fi = self.constant.imag()
                self.constant = Fr - 1j * Fi

        if (('*' in self.affstring()) or ('/' in self.affstring())
                or ('+' in self.affstring()) or ('-' in self.affstring())):
            self.string = '( ' + self.string + ' ).conj'
        else:
            self.string += '.conj' 

def eval(self, ind=None):
        if self.constant is None:
            val = spmatrix([], [], [], (self.size[0] * self.size[1], 1))
        else:
            val = self.constant
        if self.is0():
            return cvx.matrix(val, self.size)

        for k in self.factors:
            # ignore this factor if the coef is 0
            if not(self.factors[k]):
                continue
            if ind is None:
                if not k.value is None:
                    if k.vtype == 'symmetric':
                        val = val + self.factors[k] * svec(k.value)
                    else:
                        val = val + self.factors[k] * k.value[:]
                else:
                    raise Exception(k + ' is not valued')
            else:
                if ind in k.value_alt:
                    if k.vtype == 'symmetric':
                        val = val + self.factors[k] * svec(k.value_alt[ind])
                    else:
                        val = val + self.factors[k] * k.value_alt[ind][:]
                else:
                    raise Exception(
                        k + ' does not have a value for the index ' + str(ind))
        return cvx.matrix(val, self.size) 

def __init__(self, factors=None, constant=None,
                 size=(1, 1),
                 string='0'
                 ):
        if factors is None:
            factors = {}
        Expression.__init__(self, string)
        self.factors = factors
        """
                dictionary storing the matrix of coefficients of the linear
                part of the affine expressions. The matrices of coefficients
                are always stored with respect to vectorized variables (for the
                case of matrix variables), and are indexed by instances
                of the class :class:`Variable<picos.Variable>`.
                """
        self.constant = constant
        """constant of the affine expression,
                stored as a :func:`cvxopt sparse matrix <cvxopt:cvxopt.spmatrix>`.
                """
        
        assert len(size)==2 and _is_integer(size[0]) and _is_integer(size[1])
        #make sure they are of class `int`, otherwise compatibility problem with py3...
        self._size = (int(size[0]), int(size[1]))
        
        """size of the affine expression"""
        # self.string=string 

def _cplx_mat_to_real_mat(M):
    """
    if M = A +iB,
    return the block matrix [A,-B;B,A]
    """
    if not(isinstance(M, cvx.base.spmatrix) or isinstance(M, cvx.base.matrix)):
        raise NameError('unexpected matrix type')
    if M.typecode == 'z':
        A = M.real()
        B = M.imag()
    else:
        A = M
        B = spmatrix([], [], [], A.size)
    return cvx.sparse([[A, B], [-B, A]]) 

def svecm1(vec, triu=False):
    if vec.size[1] > 1:
        raise ValueError('should be a column vector')
    v = vec.size[0]
    n = int(np.sqrt(1 + 8 * v) - 1) // 2
    if n * (n + 1) // 2 != v:
        raise ValueError('vec should be of dimension n(n+1)/2')
    if not isinstance(vec, cvx.spmatrix):
        vec = cvx.sparse(vec)
    I = []
    J = []
    V = []
    for i, v in zip(vec.I, vec.V):
        c = int(np.sqrt(1 + 8 * i) - 1) // 2
        r = i - c * (c + 1) // 2
        I.append(r)
        J.append(c)
        if r == c:
            V.append(v)
        else:
            if triu:
                V.append(v / np.sqrt(2))
            else:
                I.append(c)
                J.append(r)
                V.extend([v / np.sqrt(2)] * 2)
    return spmatrix(V, I, J, (n, n)) 

def diag_vect(exp):
    """
    Returns the vector with the diagonal elements of the matrix expression ``exp``

    **Example**

    >>> import picos as pic
    >>> prob=pic.Problem()
    >>> X=prob.add_variable('X',(3,3))
    >>> pic.tools.diag_vect(X)
    # (3 x 1)-affine expression: diag(X) #

    """
    from .expression import AffinExp
    (n, m) = exp.size
    n = min(n, m)
    idx = cvx.spdiag([1.] * n)[:].I
    expcopy = AffinExp(exp.factors.copy(), exp.constant, exp.size,
                       exp.string)
    proj = spmatrix([1.] * n, range(n), idx,
                        (n, exp.size[0] * exp.size[1]))
    for k in exp.factors.keys():
        expcopy.factors[k] = proj * expcopy.factors[k]
    if not exp.constant is None:
        expcopy.constant = proj * expcopy.constant
    expcopy._size = (n, 1)
    expcopy.string = 'diag(' + exp.string + ')'
    return expcopy 

def find_nearest_valid_distribution(u_alpha, kernel, initial=None, reg=0):
    """ (solution,distance_sqd)=find_nearest_valid_distribution(u_alpha,kernel):
    Given a n-vector u_alpha summing to 1, with negative terms, 
    finds the distance (squared) to the nearest n-vector summing to 1, 
    with non-neg terms. Distance calculated using nxn matrix kernel. 
    Regularization parameter reg -- 

    min_v (u_alpha - v)^\top K (u_alpha - v) + reg* v^\top v"""

    P = matrix(2 * kernel)
    n = kernel.shape[0]
    q = matrix(np.dot(-2 * kernel, u_alpha))
    A = matrix(np.ones((1, n)))
    b = matrix(1.)
    G = spmatrix(-1., range(n), range(n))
    h = matrix(np.zeros(n))
    dims = {'l': n, 'q': [], 's': []}
    solvers.options['show_progress'] = False
    solution = solvers.coneqp(
        P,
        q,
        G,
        h,
        dims,
        A,
        b,
        initvals=initial
        )
    distance_sqd = solution['primal objective'] + np.dot(u_alpha.T,
            np.dot(kernel, u_alpha))[0, 0]
    return (solution, distance_sqd) 

def scipy_sparse_to_spmatrix(A):
    """Efficient conversion from scipy sparse matrix to cvxopt sparse matrix"""
    coo = A.tocoo()
    SP = spmatrix(coo.data.tolist(), coo.row.tolist(), coo.col.tolist(), size=A.shape)
    return SP 

def compute_displacements(self, xPhys):
        # Setup and solve FE problem
        sK = ((self.KE.flatten()[np.newaxis]).T * (
            self.Emin + (xPhys)**self.penal *
            (self.Emax - self.Emin))).flatten(order='F')
        K = scipy.sparse.coo_matrix((sK, (self.iK, self.jK)),
            shape=(self.ndof, self.ndof)).tocsc()
        # Remove constrained dofs from matrix and convert to coo
        K = deleterowcol(K, self.fixed, self.fixed).tocoo()
        # Solve system
        K1 = cvxopt.spmatrix(K.data, K.row.astype(np.int), K.col.astype(np.int))
        B = cvxopt.matrix(self.f[self.free, :])
        cvxopt.cholmod.linsolve(K1, B)
        self.u[self.free, :] = np.array(B)[:, :] 

def _cplx_vecmat_to_real_vecmat(M, sym=True, times_i=False):
    """
    if the columns of M are vectorizations of matrices of the form A +iB:
    * if times_i is False (default), return vectorizations of the block matrix [A,-B;B,A]
      otherwise, return vectorizations of the block matrix [-B,-A;A,-B]
    * if sym=True, returns the columns with respect to the sym-vectorization of the variables of the LMI
    """
    if not(isinstance(M, cvx.base.spmatrix) or isinstance(M, cvx.base.matrix)):
        raise NameError('unexpected matrix type')

    if times_i:
        M = M * 1j

    mm = M.size[0]
    m = mm**0.5
    if int(m) != m:
        raise NameError('first dimension must be a perfect square')
    m = int(m)

    vv = []
    if sym:
        nn = M.size[1]
        n = nn**0.5
        if int(n) != n:
            raise NameError('2d dimension must be a perfect square')
        n = int(n)

        for k in range(n * (n + 1) // 2):
            j = int(np.sqrt(1 + 8 * k) - 1) // 2
            i = k - j * (j + 1) // 2
            if i == j:
                v = M[:, n * i + i]
            else:
                i1 = n * i + j
                i2 = n * j + i
                v = (M[:, i1] + M[:, i2]) * (1. / (2**0.5))
            vvv = _cplx_mat_to_real_mat(cvx.matrix(v, (m, m)))[:]
            vv.append([vvv])

    else:
        for i in range(M.size[1]):
            v = M[:, i]
            A = cvx.matrix(v, (m, m))
            vvv = _cplx_mat_to_real_mat(A)[:]  # TODO 0.5*(A+A.H) instead ?
            vv.append([vvv])

    return cvx.sparse(vv) 

def _svecm1_identity(vtype, size):
    """
    row wise svec-1 transformation of the
    identity matrix of size size[0]*size[1]
    """
    if vtype in ('symmetric',):
        s0 = size[0]
        if size[1] != s0:
            raise ValueError('should be square')
        I = range(s0 * s0)
        J = []
        V = []
        for i in I:
            rc = (i % s0, i // s0)
            (r, c) = (min(rc), max(rc))
            j = c * (c + 1) // 2 + r
            J.append(j)
            if r == c:
                V.append(1)
            else:
                V.append(1 / np.sqrt(2))
        idmat = spmatrix(V, I, J, (s0 * s0, s0 * (s0 + 1) // 2))
    elif vtype == 'antisym':
        s0 = size[0]
        if size[1] != s0:
            raise ValueError('should be square')
        I = []
        J = []
        V = []
        k = 0
        for j in range(1, s0):
            for i in range(j):
                I.append(s0 * j + i)
                J.append(k)
                V.append(1)
                I.append(s0 * i + j)
                J.append(k)
                V.append(-1)
                k += 1
        idmat = spmatrix(V, I, J, (s0 * s0, s0 * (s0 - 1) // 2))
    else:
        sp = size[0] * size[1]
        idmat = spmatrix([1] * sp, range(sp), range(sp), (sp, sp))

    return idmat 

def ltrim1(vec, uptri=True,offdiag_fact=1.):
    """
    If ``vec`` is a vector or an affine expression of size n(n+1)/2, ltrim1(vec) returns a (n,n) matrix with
    the elements of vec in the lower triangle.
    If ``uptri == False``, the upper triangle is 0, otherwise the upper triangle is the symmetric of the lower one.
    """
    if vec.size[1] > 1:
        raise ValueError('should be a column vector')
    from .expression import AffinExp
    v = vec.size[0]
    n = int(np.sqrt(1 + 8 * v) - 1) // 2
    if n * (n + 1) // 2 != v:
        raise ValueError('vec should be of dimension n(n+1)/2')
    if isinstance(vec, cvx.matrix) or isinstance(vec, cvx.spmatrix):
        if not isinstance(vec, cvx.matrix):
            vec = cvx.matrix(vec)
        M = cvx.matrix(0., (n, n))
        r = 0
        c = 0
        for v in vec:
            if r == n:
                c += 1
                r = c
            if r!=c:
                v *= offdiag_fact
            M[r, c] = v
            if r > c and uptri:
                M[c, r] = v
            r += 1

        return M
    elif isinstance(vec, AffinExp):
        I, J, V = [], [], []
        r = 0
        c = 0
        for i in range(v):
            if r == n:
                c += 1
                r = c
            I.append(r + n * c)
            J.append(i)
            V.append(1)
            if r > c and uptri:
                I.append(c + n * r)
                J.append(i)
                V.append(1)
            r += 1
        H = spmatrix(V, I, J, (n**2, v))
        Hvec = H * vec
        newfacs = Hvec.factors
        newcons = Hvec.constant
        if uptri:
            return AffinExp(newfacs, newcons, (n, n),
                            'ltrim1_sym(' + vec.string + ')')
        else:
            return AffinExp(newfacs, newcons, (n, n),
                            'ltrim1(' + vec.string + ')')
    else:
        raise Exception('expected a cvx vector or an affine expression') 

def diag(exp, dim=1):
    r"""
    if ``exp`` is an affine expression of size (n,m),
    ``diag(exp,dim)`` returns a diagonal matrix of size ``dim*n*m`` :math:`\times` ``dim*n*m``,
    with ``dim`` copies of the vectorized expression ``exp[:]`` on the diagonal.

    In particular:

      * when ``exp`` is scalar, ``diag(exp,n)`` returns a diagonal
        matrix of size :math:`n \times n`, with all diagonal elements equal to ``exp``.

      * when ``exp`` is a vector of size :math:`n`, ``diag(exp)`` returns the diagonal
        matrix of size :math:`n \times n` with the vector ``exp`` on the diagonal


    **Example**

    >>> import picos as pic
    >>> prob=pic.Problem()
    >>> x=prob.add_variable('x',1)
    >>> y=prob.add_variable('y',1)
    >>> pic.tools.diag(x-y,4)
    # (4 x 4)-affine expression: Diag(x -y) #
    >>> pic.tools.diag(x//y)
    # (2 x 2)-affine expression: Diag([x;y]) #

    """
    from .expression import AffinExp
    if not isinstance(exp, AffinExp):
        mat, name = _retrieve_matrix(exp)
        exp = AffinExp({}, constant=mat[:], size=mat.size, string=name)
    (n, m) = exp.size
    expcopy = AffinExp(exp.factors.copy(), exp.constant, exp.size,
                       exp.string)
    idx = cvx.spdiag([1.] * dim * n * m)[:].I
    for k in exp.factors.keys():
        # ensure it's sparse
        mat = cvx.sparse(expcopy.factors[k])
        I, J, V = list(mat.I), list(mat.J), list(mat.V)
        newI = []
        for d in range(dim):
            for i in I:
                newI.append(idx[i + n * m * d])
        expcopy.factors[k] = spmatrix(
            V * dim, newI, J * dim, ((dim * n * m)**2, exp.factors[k].size[1]))
    expcopy.constant = cvx.matrix(0., ((dim * n * m)**2, 1))
    if not exp.constant is None:
        for k, i in enumerate(idx):
            expcopy.constant[i] = exp.constant[k % (n * m)]
    expcopy._size = (dim * n * m, dim * n * m)
    expcopy.string = 'Diag(' + exp.string + ')'
    return expcopy 

def set_sparse_upper(self, indices, bnds):
        """
        sets the upper bound bnds[i] to the index indices[i] of the variable.
        For a symmetric matrix variable, bounds on elements in the upper triangle are ignored.

        :param indices: list of indices, given as integers (column major order) or tuples (i,j).
        :type indices: ``list``
        :param bnds: list of upper bounds.
        :type lower: ``list``

        .. Warning:: This function does not modify the existing bounds on elements other
                     than those specified in ``indices``.
        """
        if self.vtype in ('hermitian', 'complex'):
            raise Exception('lower bound not supported for complex variables')
        s0 = self.size[0]
        vv = []
        ii = []
        jj = []
        for idx, up in zip(indices, bnds):
            if isinstance(idx, int):
                idx = (idx % s0, idx // s0)
            if self.vtype in ('symmetric',):
                (i, j) = idx
                if i > j:
                    ii.append(i)
                    jj.append(j)
                    vv.append(up)
                    ii.append(j)
                    jj.append(i)
                    vv.append(up)
                elif i == j:
                    ii.append(i)
                    jj.append(i)
                    vv.append(up)
            else:
                ii.append(idx[0])
                jj.append(idx[1])
                vv.append(up)
        spUP = spmatrix(vv, ii, jj, self.size)
        if self.vtype in ('symmetric',):
            spUP = svec(spUP)
        for i, j, v in zip(spUP.I, spUP.J, spUP.V):
            ii = s0 * j + i
            bil, biu = self.bnd.get(ii, (None, None))
            self.bnd._set(ii, (bil, v))
        if ('nonpositive' in self._bndtext):
            self._bndtext.replace('nonpositive', 'bounded above')
        elif ('above' not in self._bndtext) and ('upper' not in self._bndtext):
            self._bndtext += ', some upper bounds'

        for solver in self.passed:
            texteval = 'self.parent_problem.reset_' + solver + '_instance()'
            eval(texteval) 

def cvxopt_solve_qp(P, q, G=None, h=None, A=None, b=None, solver=None,
                    initvals=None, verbose=False):
    """
    Solve a Quadratic Program defined as:

    .. math::

        \\begin{split}\\begin{array}{ll}
        \\mbox{minimize} &
            \\frac{1}{2} x^T P x + q^T x \\\\
        \\mbox{subject to}
            & G x \\leq h                \\\\
            & A x = h
        \\end{array}\\end{split}

    using `CVXOPT <http://cvxopt.org/>`_.

    Parameters
    ----------
    P : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Symmetric quadratic-cost matrix.
    q : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Quadratic-cost vector.
    G : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear inequality matrix.
    h : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear inequality vector.
    A : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear equality constraint matrix.
    b : numpy.array, cvxopt.matrix or cvxopt.spmatrix
        Linear equality constraint vector.
    solver : string, optional
        Set to 'mosek' to run MOSEK rather than CVXOPT.
    initvals : numpy.array, optional
        Warm-start guess vector.
    verbose : bool, optional
        Set to `True` to print out extra information.

    Returns
    -------
    x : array, shape=(n,)
        Solution to the QP, if found, otherwise ``None``.

    Note
    ----
    CVXOPT only considers the lower entries of `P`, therefore it will use a
    wrong cost function if a non-symmetric matrix is provided.
    """
    options['show_progress'] = verbose
    args = [cvxopt_matrix(P), cvxopt_matrix(q)]
    kwargs = {'G': None, 'h': None, 'A': None, 'b': None}
    if G is not None:
        kwargs['G'] = cvxopt_matrix(G)
        kwargs['h'] = cvxopt_matrix(h)
    if A is not None:
        kwargs['A'] = cvxopt_matrix(A)
        kwargs['b'] = cvxopt_matrix(b)
    sol = qp(*args, solver=solver, initvals=initvals, **kwargs)
    if 'optimal' not in sol['status']:
        return None
    return array(sol['x']).reshape((q.shape[0],))