Python cvxopt.sparse() Examples

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 __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 _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 svec(mat, ignore_sym=False):
    """
    returns the svec representation of the cvx matrix ``mat``.
    (see `Dattorro, ch.2.2.2.1 <http://meboo.convexoptimization.com/Meboo.html>`_)

    If ``ignore_sym = False`` (default), the function raises an Exception if ``mat`` is not symmetric.
    Otherwise, elements in the lower triangle of ``mat`` are simply ignored.
    """
    if not isinstance(mat, cvx.spmatrix):
        mat = cvx.sparse(mat)

    s0 = mat.size[0]
    if s0 != mat.size[1]:
        raise ValueError('mat must be square')

    I = []
    J = []
    V = []
    for (i, j, v) in zip((mat.I), (mat.J), (mat.V)):
        if not ignore_sym:
            if abs(mat[j, i] - v) > 1e-6:
                raise ValueError('mat must be symmetric')
        if i <= j:
            isvec = j * (j + 1) // 2 + i
            J.append(0)
            I.append(isvec)
            if i == j:
                V.append(v)
            else:
                V.append(np.sqrt(2) * v)

    return spmatrix(V, I, J, (s0 * (s0 + 1) // 2, 1)) 

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 _utri(mat):
    """
    return elements of the (strict) upper triangular part of a cvxopt matrix
    """
    m, n = mat.size
    if m != n:
        raise ValueError('mat must be square')
    v = []
    for j in range(1, n):
        for i in range(j):
            v.append(mat[i, j])
    return cvx.sparse(v) 

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 __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 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 _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 __xor__(self, fact):
        """hadamard (elementwise) product"""
        selfcopy = self.copy()
        if isinstance(fact, AffinExp):
            if fact.isconstant():
                fac, facString = cvx.sparse(fact.eval()), fact.string
            else:
                if self.isconstant():
                    return fact ^ self
                else:
                    raise Exception('not implemented')
        else:
            fac, facString = _retrieve_matrix(fact, self.size[0])
        if fac.size == (1, 1) and selfcopy.size[0] != 1:
            fac = fac[0] * cvx.spdiag([1.] * selfcopy.size[0])
        if self.size == (1, 1) and fac.size[1] != 1:
            oldstring = selfcopy.string
            selfcopy = selfcopy.diag(fac.size[1])
            selfcopy.string = oldstring
        if selfcopy.size[0] != fac.size[0] or selfcopy.size[1] != fac.size[1]:
            raise Exception('incompatible dimensions')
        mm, nn = selfcopy.size
        bfac = spmatrix([], [], [], (mm * nn, mm * nn))
        for i, j, v in zip(fac.I, fac.J, fac.V):
            bfac[j * mm + i, j * mm + i] = v
        for k in selfcopy.factors:
            newfac = bfac * selfcopy.factors[k]
            selfcopy.factors[k] = newfac
        if selfcopy.constant is None:
            newfac = None
        else:
            newfac = bfac * selfcopy.constant
        selfcopy.constant = newfac
        """
                #the following removes 'I' from the string when a matrix is multiplied
                #by the identity. We leave the 'I' when the factor of identity is a scalar
                if len(facString)>0:
                        if facString[-1]=='I' and (len(facString)==1
                                 or facString[-2].isdigit() or facString[-2]=='.') and (
                                 self.size != (1,1)):
                                facString=facString[:-1]
		"""
        sstring = selfcopy.affstring()
        if len(facString) > 0:
            if ('+' in sstring) or ('-' in sstring):
                sstring = '( ' + sstring + ' )'
            if ('+' in facString) or ('-' in facString):
                facString = '( ' + facString + ' )'

            selfcopy.string = facString + '∘' + sstring

        return selfcopy 

def __rshift__(self, exp):

        if isinstance(exp, AffinExp):
            n = exp.size[0] * exp.size[1]
            if self.truncated:
                if self.nonneg:
                    if self.radius <= 1:
                        aff = (-exp[:]) // (1 | exp[:])
                        rhs = cvx.sparse([0] * n + [self.radius])
                        if self.radius == 1:
                            simptext = ' in standard simplex'
                        else:
                            simptext = ' in simplex of radius ' + \
                                str(self.radius)
                    else:
                        aff = (exp[:]) // (-exp[:]) // (1 | exp[:])
                        rhs = cvx.sparse([1] * n + [0] * n + [self.radius])
                        simptext = ' in truncated simplex of radius ' + \
                            str(self.radius)
                    cons = (aff <= rhs)
                    cons.myconstring = exp.string + simptext
                else:
                    from .problem import Problem
                    Ptmp = Problem()
                    v = Ptmp.add_variable('v', n)
                    Ptmp.add_constraint(exp[:] < v)
                    Ptmp.add_constraint(-exp[:] < v)
                    Ptmp.add_constraint((1 | v) < self.radius)
                    if self.radius > 1:
                        Ptmp.add_constraint(v < 1)
                    constring = '||' + exp.string + \
                        '||_{infty;1} <= {1;' + str(self.radius) + '}'
                    cons = Sym_Trunc_Simplex_Constraint(
                        exp, self.radius, Ptmp, constring)
            else:
                if self.nonneg:
                    aff = (-exp[:]) // (1 | exp[:])
                    rhs = cvx.sparse([0] * n + [self.radius])
                    cons = (aff <= rhs)
                    if self.radius == 1:
                        cons.myconstring = exp.string + ' in standard simplex'
                    else:
                        cons.myconstring = exp.string + \
                            ' in simplex of radius ' + str(self.radius)
                else:
                    cons = norm(exp, 1) < self.radius
            return cons
        else:  # constant
            term, termString = _retrieve_matrix(exp, None)
            exp2 = AffinExp(
                factors={},
                constant=term[:],
                size=term.size,
                string=termString)
            return self >> exp2