base.py

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

"""CVXMOD version 0.4.5, 2008-06-23."""

# Copyright (C) 2006-2008 Jacob Mattingley and Stephen Boyd.
#
# This file is part of CVXMOD.
#
# CVXMOD is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
#
# CVXMOD is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program. If not, see <http://www.gnu.org/licenses/>.

# jem
# dev at jaboc net
# jem appearing at any point is a marker that the code needs inspection.

from __future__ import division

import __builtin__
import cvxopt.solvers
import cvxopt.lapack
import cvxopt.base
from cvxopt.base import matrix, spmatrix, sparse

from math import log10, ceil
from cvxmod.errors import *

import cvxopt.info
import sys
import copy

inf = float(1e300000) # handle multiple platforms.
eps = 1e-5 # jem.

OPTVAR = 'optvar'
PARAM = 'param'
DIM = 'dim'


# jem. split this into stuff that should be imported by atoms and stuff that
# should be imported in general.
__all__ = ["affine", "classify", "cols", "concathoriz", "concatvert",
           "concave", "convex", "cvx", "decreasing", "diag", "dim", "dummyvar",
           "elementwise", "etranspose", "etp", "expandtr", "eye", "function",
           "getassertions", "getdims", "getoptvars", "getparams",
           "getquadmult", "getstdforms", "getdimmult", "getvarmult",
           "getvecmult", "equiv", "increasing", "unvec", "is1x1", "isaffine",
           "isconcave", "isconvex", "isdecreasing", "isincreasing", "isneg",
           "isnsd", "isoptvar", "isparam", "ispos", "ispsd", "issymm",
           "iszero", "lhexp", "matrix", "matrixfunction", "maximize",
           "minimize", "multfunction", "multiarg", "negative", "nnz", "nnzpc",
           "norm2", "ones", "optvar", "OPTVAR", "param", "PARAM", "positive",
           "printval", "problem", "replacevars", "rhexp", "rows", "size",
           "sparse", "speye", "spmatrix", "stdformlin", "stdformsoc",
           "stdstruct", "stroptvars", "symbol", "trace", "transpose", "tp",
           "trobj", "unitvec", "value", "vec", "vectorize", "version", "zeros",
           "zm", "issparse", "nzentries", "isnonzero"]
# remove getvarmult, getquadmult, info later. change name of seed.
# remove stdformsoc. remove stdstruct. remove trobj.

def version():
    print __doc__
    print 'Using CVXOPT version', cvxopt.info.version + '.'
    print 'Platform is ' + sys.platform + '.'
    print 'Using Python version %s, located at %s.' % (sys.version[:5], sys.executable)

class lhsrhs(object):
    def getoptvars(self):
        return getoptvars(self.lhs) | getoptvars(self.rhs)

    def getassertions(self):
        a = getassertions(self.lhs) + getassertions(self.rhs)
        # jem add alternate assertions: if either is 1x1, fine; otherwise,
        # match up sizes.
        return a

    def getstdforms(self):
        return getstdforms(self.lhs) | getstdforms(self.rhs)

    def getparams(self):
        return getparams(self.lhs) | getparams(self.rhs)

    def getdims(self):
        return getdims(self.lhs) | getdims(self.rhs)

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

        if isoptvar(self.rhs) and self.rhs in set(d.keys()):
            self.rhs = d[self.rhs]
        else:
            replacevars(self.rhs, d)
    def _getconvex(self):
        return isconvex(self.lhs) and isconvex(self.rhs)
    convex = property(_getconvex)

    def _getconcave(self):
        return isconcave(self.lhs) and isconcave(self.rhs)
    concave = property(_getconcave)

    def _getdecreasing(self):
        return isdecreasing(self.lhs) and isdecreasing(self.rhs)
    decreasing = property(_getdecreasing)

    def _getincreasing(self):
        return isincreasing(self.lhs) and isincreasing(self.rhs)
    increasing = property(_getincreasing)
            
    def _getpos(self):
        return ispos(self.lhs) and ispos(self.rhs)
    pos = property(_getpos)
            
    def _getneg(self):
        return isneg(self.lhs) and isneg(self.rhs)
    neg = property(_getneg)

class convex(object):
    """General holder class for convex functions."""
    convfn = True

class concave(object):
    """General holder class for concave functions."""
    concfn = True

class affine(convex, concave):
    """General holder class for affine functions."""
    # jem. need to update with (possibly) mapping functions or something.
    convfn = True
    concfn = True

class increasing(object):
    """General holder class for increasing functions."""
    incfn = True

class decreasing(object):
    """General holder class for decreasing functions."""
    decfn = True

class positive(object):
    posfn = True

class negative(object):
    negfn = True

class psd(object):
    psd = True
    nsd = False
    symm = True

class nsd(object):
    psd = False
    nsd = True
    symm = True

class symbol(affine):
    """A general holder for optvars or params involved in optimization."""
    increasing = True
    convex = True
    concave = True
    def __init__(self, name=None, rs=None, cs=None, value=None, role=PARAM,
                 pos=False, neg=False, symm=False, psd=False,
                 nsd=False, lower=None, upper=None):
        # Prevent some of these attributes from changing?
        # jem.
        # optional warning levels?
        #if not isinstance(name, str):
        #    Jwarn("using string '%s' as name" % str(name))
        if name is None:
            self.name = 'noname' + str(nonamecount())
        else:
            self.name = name

        if rs is None and cs is None and value is not None:
            rs = rows(value)
            cs = cols(value)
        else:
            if rs is None:
                rs = 1
            if cs is None:
                cs = 1

        if not is1x1(rs):
            raise TypeError("rs must be a 1x1 quantity")
        self.rows = rs

        if not is1x1(rs):
            raise TypeError("cs must be a 1x1 quantity")
        self.cols = cs

        # jem: later change these warnings to errors?
        if getparams(rs):
            Jwarn("rows should not contain parameters: use dims instead (got %s)" % str(rs))

        if getparams(cs):
            Jwarn("cols should not contain parameters: use dims instead (got %s)" % str(cs))

        self.role = role

        self.pos = pos
        self.neg = neg
        self.symm = symm
        self.psd = psd
        self.nsd = nsd
        self.lower = lower
        self.upper = upper
        
        # jem. add assertions about size. deal with parametric size.
        if value is not None:
            if isdim(rs):
                # jem: add assertion here.
                pass
            #elif rs is not rows(value):
            #    raise ValueError("value must match in size")

            if isdim(cs):
                # jem: add assertion here.
                pass
            #elif cs is not cols(value):
            #    raise ValueError("value must match in size")

        if value is not None:
            self.value = value__(value)
        else:
            self.value = None

        self._vecvar = None

    def __repr__(self):
        if isdim(self):
            return "<dim symbol %s>" % self.name

        v = self.getattrs()
        if v:
            extras = '(' + ', '.join(v) + ') '
        else:
            extras = ''

        return "<%sx%s %s%s symbol %s>" % (withbrackets(self.rows),
                                            withbrackets(self.cols), extras,
                                            self.role, self.name)

    def __getitem__(self, k):
        return symbslice(self, k)

    def getvarmult(self, var):
        if self.role in [OPTVAR, 'unitvec']:
            if self is var:
                return eye(rows(var), True)
            else:
                return 0
        elif isparam(self) and var is None:
            return self
        else:
            return 0

    def getdimmult(self, dim):
        if isdim(self) and self is dim:
            return 1
        else:
            return 0

    def getvecmult(self, var):
        if self.role in [OPTVAR, 'unitvec']:
            if self is var:
                return eye(rows(var)*cols(var), True)
            else:
                return 0
        elif isparam(self) and var is None:
            return vec(self)
        else:
            return 0

    def __mul__(self, other):
        return multiply(self, other)
    
    def __rmul__(self, other):
        return multiply(other, self)

    def __rdiv__(self, other):
        return other*(self**-1)

    def __add__(self, other):
        return addup(self, other)

    def __radd__(self, other):
        return addup(self, other)

    def __sub__(self, other):
        return addup(self, -other)
    
    def __rsub__(self, other):
        return addup(other, -self)

    def __neg__(self):
        return multiply(-1, self)

    def __pos__(self):
        return self

    def __lt__(self, other):
        return relation(self, other, '<')

    def __le__(self, other):
        return relation(self, other, '<=')

    def __gt__(self, other):
        return relation(self, other, '>')

    def __ge__(self, other):
        return relation(self, other, '>=')

    def __eq__(self, other):
        return relation(self, other, '==')

    def __float__(self):
        if is1x1(self):
            if self.value is None:
                raise OptvarValueError
            else:
                return float(self.value)
        else:
            raise ValueError('cannot convert %s to float' % self.name)

    def __int__(self):
        if is1x1(self):
            if self.value is None:
                raise OptvarValueError
            else:
                return int(self.value)
        else:
            raise ValueError('cannot convert %s to int' % self.name)

    def __str__(self):
        return self.name

    def __pow__(self, other):
        if other is -1:
            return inv(self)

    def getattrs(self):
        v = []
        if self.symm:
            v.append('symmetric')
        if self.psd:
            v.append('psd')
        if self.nsd:
            v.append('nsd')
        if self.neg:
            v.append('neg')
        if self.pos:
            v.append('pos')

        return v

    def getoptvars(self):
        if isoptvar(self):
            return set((self,))
        else:
            return set()

    def getassertions(self):
        a = [ensureint(rows(self)), ensureint(cols(self))]
        if self.symm:
            a.append(rows(self) == cols(self))

        return a

    def getparams(self):
        if isparam(self):
            return set((self,)) | getparams(self.rows) | getparams(self.cols)
        else:
            return set()

    def getdims(self):
        if isdim(self):
            return set((self,))
        else:
            return getdims(self.rows) | getdims(self.cols)

    def cvxdeclare(self):
        if isoptvar(self):
            s = ""
            if is1x1(self):
                s += "variable %s" % str(self)
            elif self.cols is 1:
                s += "variable %s(%s)" % (str(self), str(self.rows))
            else:
                s += "variable %s(%s,%s)" % (str(self), str(self.rows),
                                               str(self.cols))
            if self.symm:
                s += " symmetric"
            if self.psd:
                s += "\n    %s == semidefinite(%s)" % (str(self), str(self.rows))
            return s
        else:
            raise OptimizationError('cannot call cvxdeclare on non-optvar')

    def matlab(self):
        if isparam(self):
            s = self.name + ' = '
            if isinstance(self.value, (matrix, spmatrix)):
                s += '['
                offset = ' '*len(s)
                for i in xrange(value(rows(self))):
                    for j in xrange(value(cols(self))):
                        s += '% .5f'.rjust(6) % self.value[i,j]
                    s += ';\n'
                    s += offset

                s = s[:-(2+len(offset))] + '];\n'
            else:
                s += '% .4g'.rjust(5)
                s += ';\n' % self.value[i,j]

            return s
        else:
            return OptimizationError('cannot call matlab on non-param')

    def _getdecreasing(self):
        return not isoptvar(self)