base.py

CVXMOD is a Python-based tool for expressing and solving convex optimization problems.

        return expandtr(self.lhs)*expandtr(self.rhs)

    def getquadmult(self, var):
        if var not in getoptvars(self):
            return 0

        c = self

        lhs = 1
        rhs = 1
        while isinstance(c, multfunction):
            if var in getoptvars(c.lhs) and var not in getoptvars(c.rhs):
                rhs = c.rhs * rhs
                c = c.lhs
            elif var not in getoptvars(c.lhs) and var in getoptvars(c.rhs):
                lhs = lhs * c.lhs
                c = c.rhs
            elif var in getoptvars(c.lhs) and var in getoptvars(c.rhs):
                lhs, rhs = lhs*c.lhs, c.rhs*rhs

                lhm = getvarmult(etranspose(lhs), var)
                rhm = getvarmult(rhs, var)

                return etranspose(lhm)*rhm
            else:
                raise OptimizationError("shouldn't be here")

        return 0

    def vectorize(self):
        V = getoptvars(self.lhs)
        if len(V) > 1:
            raise OptimizationError('cannot vectorize: too many variables in'
                                    'expression')
        n = value(cols(self))
        lhs = self.lhs
        if V:
            v = V.pop()
            if not is1x1(v) or getoptvars(self.rhs):
                raise OptimizationError('cannot vectorize when non-1x1'
                                        'variable is on lhs')

        return [self.lhs*x for x in vectorize(self.rhs)]
    def cvx(self):
        if hasattr(self.lhs, 'brackets') and self.lhs.brackets:
            l = '(%s)' % cvx(self.lhs)
        else:
            l = cvx(self.lhs)

        if hasattr(self.rhs, 'brackets') and self.rhs.brackets:
            r = '(%s)' % cvx(self.rhs)
        else:
            r = cvx(self.rhs)
            
        return l + '*' + r

    def _getsymm(self):
        return is1x1(self.lhs) and issymm(self.rhs)
    symm = property(_getsymm)

    def nnz(self):
        if is1x1(self.lhs):
            return nnz(self.rhs)
        elif is1x1(self.rhs):
            return nnz(self.lhs)
        else:
            return rows(self)*cols(self)

    def nzentries(self):
        if is1x1(self.lhs):
            d = nzentries(self.rhs)
            for k in d.keys():
                d[k] = self.lhs*d[k]

            return d
        else:
            raise NotImplementedError

def multiply(lhs, rhs):
    # jemjemjem size checking should be in this function, not in
    # multfunction!!!

    # jem: speed optimized. Be gentle.

    if iszero(lhs):
        return 0
    elif iszero(rhs):
        return 0
    elif lhs is 1 or lhs is 1.0:
        return rhs
    elif rhs is 1 or rhs is 1.0:
        return lhs
    elif not is1x1(lhs) and isinstance(rhs, eyesymbol):
        return lhs # jemjem. need size checking here.
    elif not is1x1(rhs) and isinstance(lhs, eyesymbol):
        return rhs # jemjem. need size checking here.
    else:
        a = multfunction(lhs, rhs)
        # check that we don't have --x == 1*x.
        if is1x1(a.lhs) and (a.lhs is 1 or a.lhs is 1.0):
            return a.rhs
        else:
            return a
    
def addup(lhs, rhs):
    if iszero(lhs):
        return rhs
    elif iszero(rhs):
        return lhs
    else:
        return addfunction(lhs, rhs)

def JGen(i = 0):
    """Returns a generator of incrementing integers."""
    while 1:
        yield i
        i += 1

def tcount(init = JGen(1)):
    return init.next()

def nonamecount(init = JGen(1)):
    return init.next()

class norm2function(function, convex, positive):
    name = 'norm2'
    def __init__(self, arg):
        self.arg = arg
        self.rows = 1
        self.cols = 1

    def _getvalue(self):
        return sum(value(self.arg)**2)**0.5
    value = property(_getvalue)

    def _getincfn(obj):
        return ispos(obj.arg)
    incfn = property(_getincfn)

    def _getdecfn(obj):
        return isneg(obj.arg)
    decfn = property(_getdecfn)

    def cvx(self):
        return 'norm(%s, 2)' % str(self.arg)

    def epiorhypo(self):
        t = dummyvar(1, 1)
        t2 = dummyvar(size(self.arg))
        if isoptvar(self.arg):
            s = (norm2(self.arg) <= t).split()
        elif isaffine(self.arg):
            s = (self.arg == t2).split()
            s += (norm2(t2) <= t).split()
        else:
            if self.incfn:
                s = (self.arg <= t2).split()
            elif self.decfn:
                s = (self.arg >= t2).split()
            else:
                raise ConvexityError('invalid norm2 composition')
            s += (norm2(t2) <= t).split()
        return (t, s)

def norm2(obj):
    if isinstance(obj, (int, float, matrix, spmatrix)):
        return sum(obj**2)**0.5
    else:
        # jem. this (and possibly other situations) will be problematic with
        # params. need some kind of conditional comparison / put on stack
        # for later assertion.
        if cols(obj) is not 1:
            raise DimError('norm2 arg must be a column vector')
        else:
            return norm2function(obj)

class relation(lhsrhs):
    # split this according to sign [possibly w/ superclass]? then don't have to
    # use self.sign as an attribute.
    # should this have rows and columns???
    def __init__(self, lhs, rhs, sign):
        # If either is a scalar, no worries. If they're both non-1x1, they
        # had better be the same (at this point, anyway).
        if is1x1(lhs) or is1x1(rhs) or equivdims(rows(lhs), rows(rhs)) and \
           equivdims(cols(lhs), cols(rhs)):
            self.lhs = lhs
            self.rhs = rhs
            self.sign = sign
            if is1x1(lhs):
                self.rows = rows(rhs)
                self.cols = cols(rhs)
            else:
                self.rows = rows(lhs)
                self.cols = cols(lhs)
        else:
            raise DimError

    def _getvalue(self):
        if self.sign == '==':
            # could have more sophisticated checking in here.
            a = value(self.lhs - self.rhs)
            if not isinstance(a, (matrix, spmatrix)):
                a = matrix(a, tc='d')
            return cvxopt.blas.nrm2(a) <= (eps**2)*rows(a) + 1e-7
        else:
            # jem. deal with bug of no proper comparison for matrices.
            a = self.stdlhs()
            A = value(a)
            if isinstance(A, (int, float)):
                return A < eps
            for i in xrange(len(A)):
                if A[i] > eps:
                    return False
            return True
    value = property(_getvalue)

    def __repr__(self):
        if self.sign == '==':
            return "<%sx%s equality %s %s %s; optvars: %s>" % \
                (str(rows(self)), str(cols(self)), str(self.lhs), self.sign,
                str(self.rhs), stroptvars(self))
        else:
            return "<%sx%s inequality %s %s %s; optvars: %s>" % \
                (str(rows(self)), str(cols(self)), str(self.lhs), self.sign,
                str(self.rhs), stroptvars(self))

    def __str__(self):
        return "%s %s %s" % (str(self.lhs), self.sign, str(self.rhs))

    def cvx(self):
        return "%s %s %s" % (cvx(self.lhs), self.sign, cvx(self.rhs))

    def latex(self):
        s = self.sign
        if s == '<=':
            s = r'\leq'
        elif s == '>=':
            s = r'\geq'

        return "%s %s %s" % (latex(self.lhs), s, latex(self.rhs))

    def _getaffine(self):
        return isaffine(self.lhs - self.rhs)
    affine = property(_getaffine)

    def _getconvex(self):
        # This isn't really defined?
        if self.sign == '==':
            return isaffine(self.lhs) and isaffine(self.rhs)
        elif self.sign in ('<=', '<'):
            return isconvex(self.lhs - self.rhs)
        else:
            return isconvex(self.rhs - self.lhs)
    convex = property(_getconvex)

    def __nonzero__(self):
        raise OptimizationError("you should use `is', not `==' for comparing "
              "expressions/optvars.\nalso don't use bool(a >= 0), for example")

    def stdlhs(self):
        if self.sign == '<':
            return self.lhs - self.rhs
        elif self.sign == '<=':
            return self.lhs - self.rhs
        elif self.sign == '>':
            return self.rhs - self.lhs
        elif self.sign == '>=':
            return self.rhs - self.lhs
        else:
            return self.lhs - self.rhs

    def split(self):
        if isaffine(self):
            return stdstruct(lin=[self])
        elif isinstance(self.lhs, norm2function) and isoptvar(self.lhs.arg) and isoptvar(self.rhs):
            return stdstruct(socone=[self])
        else:
            (t, s) = self.stdlhs().epiorhypo()
            s.lin.append(t <= 0)
            return s

class problem(object):
    """problem([objective[, constr[, assertions]]]) -> convex optimization
    problem"""
    # jem. add affine property for problem.
    def __init__(self, objective=None, constr=None, assertions=None):
        self.objective = objective
        if constr is None:
            self.constr = []
        else:
            self.constr = constr
        if assertions is None:
            self.assertions = []
        else:
            self.assertions = assertions

    def __str__(self):
        if self.constr:
            s = "%s,\nsubject to:\n" % self.strobjective()
            for x in self.constr:
                s += '  ' + str(x) + '\n'
            s = s[:-1]
        elif self.objective is not None:
            s = self.strobjective()
        else:
            return 'empty problem.'

        s += '\n\noptimization variables:\n'
        for x in bylowerstr(getoptvars(self)):
            s += '  %s\n' % repr(x)
        s = s[:-1]

        if getparams(self):
            s += '\nparameters:\n'
            for x in bylowerstr(getparams(self)):
                s += '  %s\n' % repr(x)
            s = s[:-1]

        if getdims(self):
            s += '\ndimensions:\n'
            for x in bylowerstr(getdims(self)):
                s += '  %s\n' % repr(x)
            s = s[:-1]

        return s

    def __repr__(self):
        if len(self.constr) == 0:
            constr = 'no contraints'
        elif len(self.constr) == 1:
            constr = '1 constraint'
        else:
            constr = '%d constraints' % len(self.constr)

        if len(self.assertions) == 0:
            assertions = 'no assertions'
        elif len(self.assertions) == 1:
            assertions = '1 assertion'
        else:
            assertions = '%d assertions' % len(self.assertions)

        return "<problem with %s and %s; optvars: %s>" % \
                (constr, assertions, stroptvars(self))

    def gettype(self):
        if self.objective is None:
            return "feasibility"
        else:
            return self.objective.gettype()

    def cvx(self):
        s = ''
        s += '% CVX representation, generated by CVXMOD.\n'
        s += '% WARNING: not all atoms are supported in CVX.\n'
        s += '% WARNING: code may require modification.\n'

        if self.getparams():
            s += '\n% Parameter definitions.\n'

        for x in bylowerstr(self.getparams()):
            if x.value:
                s += x.matlab()
            else:
                pass

        if self.getdims():
            s += '\n% Dimension definitions.\n'

        for x in bylowerstr(self.getdims()):
            if x.value:
                s += x.matlab()
            else:
                pass

        s += "\ncvx_begin\n"
        for x in getoptvars(self):
            s += "    " + x.cvxdeclare() + '\n'
            if x.pos:
                s += "    " + str(x >= 0) + '\n'
            if x.neg:
                s += "    " + str(x <= 0) + '\n'

        if self.objective is not None:
            s += "    %s\n" % cvx(self.objective)
            if len(self.constr) > 0:
                s += "    subject to\n"

        for x in self.constr:
            if self.objective is None:
                # Don't indent as much.
                s += "    %s\n" % cvx(x)
            else:
                s += "        %s\n" % cvx(x)
        s += "cvx_end"
        return s

    def latex(self):
        s = ''
        s += '% Latex representation, generated by CVXMOD.\n'
        s += '% WARNING: markup may require modification.\n'

        if self.getparams():
            s += '\nLet \n'

        s += ', '.join([x.latexintroduce for x in self.getparams()]) + '.'

        s += '\nIn particular, let\n\\[\n'

        for x in self.getparams():
            if x.value:
                s += x.latexintroduce()
            else:
                pass

        s += '\nIn particular, let\n\\[\n'

        s += "\ncvx_begin\n"
        for x in self.getoptvars():
            s += "    " + x.latexdeclare() + '\n'
        if self.objective is not None:
            s += "    %s\n" % cvx(self.objective)
            if len(self.constr) > 0:
                s += "    subject to\n"

        for x in self.constr:
            if self.objective is None:
                # Don't indent.
                s += " %s\n