base.py

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

    decreasing = property(_getdecreasing)

    def _columnize(self):
        # jem verbose notes (TM).
        # could cunningly just have a columns property that automatically
        # creates these the first time, to make it transparent (and do the same
        # with transposed version, although add extra constr to self.)
        # produce (and cache) column columns to represent this matrix variable.
        # at this stage we are solving a problem, so all params should have
        # valid values.
        (m, n) = value(size(self))

        if not isoptvar(self):
            raise OptimizationError('should not be columnizing a param')
        if n == 1:
            raise OptimizationError('should not be columnizing a column')

        self._columnsymbs = []
        # we can proceed. start by making appropriate number of columns.
        # indexing from zero, because internal only.
        for i in range(n):
            self._columnsymbs.append(optvar(self.name + '_c' + str(i), m, 1))

            # this next bit is for the required extra inequalities.
            #s.lin += (X.columns[j] == \ 
            #          self.arg*unitvec(j+1, cols(self.arg))).split()

    def _getcolumns(self):
        try:
            return self._columnsymbs
        except AttributeError:
            self._columnize()
            return self._columnsymbs
    columns = property(_getcolumns)
    # jem: consider making columns, backsubcols private?

    def backsubcols(self):
        """Moves values from self.vectors to self.value."""
        (m, n) = value(size(self))
        self.value = zeros(m, n)
        for i in range(n):
            self.value[:,i] = value(self.columns[i])

    def vectorize(self):
        if isoptvar(self):
            return self.columns
        else:
            x = []
            A = value(self)
            for i in range(value(cols(self))):
                x.append(A[:,i])
            return x

    def getvecvar(self):
        # Create a vectorized version of this variable.
        if self.symm:
            # int() is to enforce data type only, in conjunction with importing
            # division from __future__.
            m = value(rows(self))
            numfreevars = int(m*(m+1)/2)
            return unsymm(value(rows(self)))*dummyvar(numfreevars)
        else:
            return dummyvar(value(rows(self)*cols(self)), 1)

    def nzentries(self):
        if isparam(self):
            d = {}
            for i in range(value(rows(self))):
                for j in range(value(cols(self))):
                    d[i,j] = self[i,j]
            return d
        else:
            raise NotImplementedError

class symbslice(symbol):
    """Sliced (indexed) symbol, e.g. x[0] or A[2,3] or X[(1,2,3),0]."""
    def __init__(self, arg, sl):
        self.arg = arg
        self.sl = sl
        # Pairmode means that sl is a sequence of indices, i.e. diag(X) would
        # have X[((0,0), (1,1), (2,2))].
        if iterable(sl) and iterable(sl[0]) and iterable(sl[1]):
            self.pairmode = True
        else:
            self.pairmode = False
        self.role = 'slice'

    def _getrows(self):
        (rowind, colind) = self.splitsl()
        if isinstance(rowind, slice):
            if colind is None:
                a = self.getslicelen(rowind, rows(self.arg)*cols(self.arg))
            else:
                a = self.getslicelen(rowind, rows(self.arg))
            if a is None:
                # jem: this is for dimensions not yet known?
                # jem: update this to be more sensible? can calculate lengths
                # in terms of the parameter (though only to some extent).
                # could add an assertion here about the eventual value of
                # symbol.
                return symbol('?')
            else:
                return a
        elif iterable(rowind):
            return len(rowind)
        else:
            return 1
    rows = property(_getrows)

    def _getcols(self):
        (rowind, colind) = self.splitsl()
        if colind is None:
            return 1
        elif self.pairmode:
            return 1
        elif isinstance(colind, slice):
            a = self.getslicelen(colind, cols(self.arg))
            if a is None:
                # jem: this is for dimensions not yet known?
                # jem: update this to be more sensible? can calculate lengths
                # in terms of the parameter (though only to some extent).
                # could add an assertion here about the eventual value of
                # symbol.
                return symbol('?')
            else:
                return a
        elif iterable(colind):
            return len(colind)
        else:
            return 1
    cols = property(_getcols)

    def getslicelen(self, s, l):
        if s.start is None:
            start = 0
        else:
            start = s.start

        if s.stop is None:
            stop = l
        elif isinstance(s.stop, int) and s.stop < 0:
            stop = l + s.stop
        else:
            stop = s.stop

        if s.step is None:
            step = 1
        elif not isinstance(s.step, int):
            return None
        elif s.step < 0:
            # Simply make positive, for the point of view of lengths.
            step = -s.step
        else:
            step = s.step

        # just give up if anything is negative.
        if (not isinstance(start, int)) or start < 0:
            raise NotImplementedError('unsupported start index')

        if step is 1:
            return stop - start
        else:
            if isinstance(start, int) and isinstance(stop, int):
                return int(ceil((stop - start) / float(step)))
            else:
                raise NotImplementedError('unsupported start or stop indices')

    def __repr__(self):
        return "<%sx%s sliced %s symbol %s>" % (withbrackets(self.rows),
                                                withbrackets(self.cols),
                                                self.arg.role, str(self))

    def splitsl(self):
        """Splits the slice into its left- and right-hand sides."""
        # If both are lists, they have to be (row, col) pairs.
        if not iterable(self.sl):
            rowind = self.sl
            colind = None
        elif iterable(self.sl) and len(self.sl) == 2 and \
                not (iterable(self.sl[0]) and iterable(self.sl[1])):
            # Not in (row, col) pairs format.
            rowind = self.sl[0]
            colind = self.sl[1]
        else:
            # Assume in (row, col) pairs format.
            (rowind, colind) = zip(*self.sl)
        return (rowind, colind)

    def __str__(self):
        if self.pairmode:
            slstr = str(self.sl)
        else:
            if iterable(self.sl):
                (rowind, colind) = self.sl
            else:
                rowind = self.sl
                colind = None

            slstr = ''
            if isinstance(rowind, slice):
                if rowind.start:
                    slstr += str(rowind.start)
                slstr += ':'
                if rowind.stop:
                    slstr += str(rowind.stop)
                if rowind.step:
                    slstr += ':' + str(rowind.step)
            else:
                slstr += str(rowind)

            if colind is not None:
                slstr += ','
                if isinstance(colind, slice):
                    if colind.start:
                        slstr += str(colind.start)
                    slstr += ':'
                    if colind.stop:
                        slstr += str(colind.stop)
                    if colind.step:
                        slstr += ':' + str(colind.step)
                else:
                    slstr += str(colind)

        return "%s[%s]" % (str(self.arg), slstr)

    def getoptvars(self):
        return getoptvars(self.arg)

    def getparams(self):
        return getparams(self.arg)

    def _getvalue(self):
        (rowind, colind) = self.splitsl()
        if colind is None:
            return value(self.arg)[rowind]
        else:
            return value(self.arg)[rowind,colind]
    value = property(_getvalue)

    def _getincreasing(self):
        return isincreasing(self.arg)
    increasing = property(_getincreasing)

    def _getdecreasing(self):
        return isdecreasing(self.arg)
    decreasing = property(_getdecreasing)

    def _getconvex(self):
        return isconvex(self.arg)
    convex = property(_getconvex)
    
    def _getconcave(self):
        return isconcave(self.arg)
    concave = property(_getconcave)

    def _getpos(self):
        return ispos(self.arg)
    pos = property(_getpos)

    def _getneg(self):
        return isneg(self.arg)
    neg = property(_getneg)

    def getvecmult(self, var):
        m = value(rows(self.arg))
        n = value(cols(self.arg))
        if var in getoptvars(self.arg):
            (rowind, colind) = self.splitsl()
            if self.pairmode:
                L = zeros(len(rowind), value(rows(self.arg)*cols(self.arg)))
                for k in range(len(rowind)):
                    L[k,rowind[k]+m*colind[k]] = 1
                return getvecmult(L*vec(self.arg), var)

            if isinstance(rowind, slice):
                rowindices = slicetoindices(rowind, rows(self.arg))
            elif iterable(rowind):
                rowindices = rowind
            else:
                rowindices = [value(rowind)]

            if colind is None:
                L = zeros(len(rowindices), value(rows(self.arg)*cols(self.arg)))
                for (i, j) in zip(range(len(rowindices)), rowindices):
                    L[i,j] = 1
                return getvecmult(L*self.arg, var)
            else:
                colindices = slicetoindices(colind, cols(self.arg))

                L = zeros(len(rowindices), value(rows(self.arg)))
                for (i, j) in zip(range(len(rowindices)), rowindices):
                    L[i,j] = 1

                R = zeros(value(cols(self.arg)), len(colindices))
                for (i, j) in zip(colindices, range(len(colindices))):
                    R[i,j] = 1

                return getvecmult(L*self.arg*R, var)
        else:
            return 0

    def getvarmult(self, var):
        return getvecmult(self, var)

    def cvx(self):
        if is1x1(self):
            (rowind, colind) = self.splitsl()
            if colind is None:
                return '%s(%d)' % (str(self.arg), rowind+1)
            else:
                return '%s(%d,%d)' % (str(self.arg), rowind+1, colind+1)
        else:
            raise NotImplementedError('cannot cvx complicated slices')

def optvar(name=None, rows=None, cols=None, value=None, role=OPTVAR,
           pos=False, neg=False, symm=False, psd=False, nsd=False, lower=None,
           upper=None):
    """Equivalent to calling symbol(), but with default role='optvar' instead
    of 'param'."""
    return symbol(name, rows, cols, value, role, pos, neg, symm, psd,
                  nsd, lower, upper)
    
def param(name=None, rows=None, cols=None, value=None, role=PARAM,
          pos=False, neg=False, symm=False, psd=False, nsd=False):
    """Equivalent to calling symbol(), with default role='param'."""
    return symbol(name, rows, cols, value, role, pos, neg, symm, psd,
                  nsd)

def dim(name=None, value=None):
    """Equivalent to calling symbol(), but with default role='dim' instead
    of 'param'."""
    # Most of the parameters are set to False---even though they are
    # technically True, they're irrelevant.
    return symbol(name, rs=1, cs=1, value=value, role='dim', pos=True,
                  neg=False, symm=False, psd=False, nsd=False)
    
class _unitvec(symbol):
    def __init__(self, i, rows):
        m = zeros(rows,1)
        m[i] = 1
        symbol.__init__(self, 'e'+str(i), rows, 1, m, 'unit vector', False, False,
                        False, False, False)
        self.i = i

    def getvarmult(self, var):
        if var is None:
            return self

def unitvec(i, rows):
    # could sneakily allow an option to have an indeterminate number of rows -
    # resolves at runtime from requirements for multiplication to be correct.
    # jem.
    if i > rows - 1:
        raise OptimizationError('unit vector index cannot exceed number'
                                       'of rows - 1')
    elif i < 0:
        raise OptimizationError('unit vector cannot be less than 0')
    else:
        if rows == 1:
            return 1
        else:
            return _unitvec(i,rows)

class objective(object):
    """Holds objectives / convex optimization goals."""
    def __repr__(self):
        return "<objective %s; optvars: %s>" % (str(self),
                                                  stroptvars(self))

    def getoptvars(self):
        return getoptvars(self.arg)

    def getassertions(self):
        return getassertions(self.arg)

    def getstdforms(self):
        return getstdforms(self.arg)
    
    def getparams(self):
        return getparams(self.arg)

    def getdims(self):
        return getdims(self.arg)

    def _getvalue(self):
        # perhaps later make this [0] to avoid returning a matrix?
        return value(self.arg)
    value = property(_getvalue)

    def classify(self):
        if isconvex(self):
            print "convex objective."
        else:
            print "non-convex objective."

    def replacevars(self, d):
        if isoptvar(self.arg) and self.arg in set(d.keys()):
            self.arg = d[self.arg]
        else:
            replacevars(self.arg, d)

class minimizeobjective(objective):
    def __init__(self, arg):
        self.arg = arg
        self.rows = rows(arg)
        self.cols = cols(arg)

    def __str__(self):
        return "minimize %s" % str(self.arg)

    def cvx(self):
        return "minimize(%s)" % cvx(self.arg)

    def latex(self):
        return r"\mbox{mimimize} & %s" % latex(self.arg)

    def gettype(self):
        return "minimization"

    def _getconvex(self):
        return isconvex(self.arg)
    convex = property(_getconvex)

    def minarg(self):
        return self.arg
        
        # jem experimental.
        t = dummyvar(size(self.arg))
        return (minimize(t), t >= self.arg)