# Source code for cvxpy.constraints.exponential

"""
Copyright 2013 Steven Diamond

This file is part of CVXPY.

CVXPY 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.

CVXPY 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 CVXPY.  If not, see <http://www.gnu.org/licenses/>.
"""

import cvxpy.settings as s
from cvxpy.error import SolverError
import cvxpy.lin_ops.lin_utils as lu
from cvxpy.lin_ops.lin_op import VARIABLE
import cvxpy.utilities.performance_utils as pu
from cvxpy.constraints.nonlinear import NonlinearConstraint
from cvxpy.constraints.utilities import format_elemwise
import numpy as np
import math

[docs]class ExpCone(NonlinearConstraint):
"""A reformulated exponential cone constraint.

Operates elementwise on :math:x, y, z.

Original cone:

.. math::

K = \\{(x,y,z) \mid y > 0, ye^{x/y} <= z\\}
\\cup \\{(x,y,z) \mid x \leq 0, y = 0, z \geq 0\\}

Reformulated cone:

.. math::

K = \\{(x,y,z) \mid y, z > 0, y\\log(y) + x \\leq y\\log(z)\\}
\\cup \\{(x,y,z) \\mid x \\leq 0, y = 0, z \\geq 0\\}

Parameters
----------
x : Variable
x in the exponential cone.
y : Variable
y in the exponential cone.
z : Variable
z in the exponential cone.
"""

def __init__(self, x, y, z, constr_id=None):
self.x = x
self.y = y
self.z = z
super(ExpCone, self).__init__(self._solver_hook,
[self.x, self.y, self.z],
constr_id)

def __str__(self):
return "ExpCone(%s, %s, %s)" % (self.x, self.y, self.z)

def __repr__(self):
return "ExpCone(%s, %s, %s)" % (self.x, self.y, self.z)

@property
def residual(self):
# TODO(akshayka): The projection should be implemented directly.
from cvxpy import Problem, Minimize, Variable, norm2, hstack
if self.x.value is None or self.y.value is None or self.z.value is None:
return None
x = Variable(self.x.shape)
y = Variable(self.y.shape)
z = Variable(self.z.shape)
constr = [ExpCone(x, y, z)]
obj = Minimize(norm2(hstack([x, y, z]) -
hstack([self.x.value, self.y.value, self.z.value])))
problem = Problem(obj, constr)
return problem.solve()

def format(self, eq_constr, leq_constr, dims, solver):
"""Formats EXP constraints for the solver.

Parameters
----------
eq_constr : list
A list of the equality constraints in the canonical problem.
leq_constr : list
A list of the inequality constraints in the canonical problem.
dims : dict
A dict with the dimensions of the conic constraints.
solver : str
The solver being called.
"""
if solver.name() == s.CVXOPT:
eq_constr += self.__CVXOPT_format[0]
elif solver.name() in [s.SCS, s.JULIA_OPT]:
leq_constr += self.__SCS_format[1]
elif solver.name() == s.ECOS:
leq_constr += self.__ECOS_format[1]
else:
raise SolverError("Solver does not support exponential cone.")
# Update dims.
dims[s.EXP_DIM] += self.num_cones()

@property
def size(self):
"""The number of entries in the combined cones.
"""
# TODO use size of dual variable(s) instead.
return sum(self.cone_sizes())

def num_cones(self):
"""The number of elementwise cones.
"""
return np.prod(self.args[0].shape, dtype=int)

def cone_sizes(self):
"""The dimensions of the exponential cones.

Returns
-------
list
A list of the sizes of the elementwise cones.
"""
return [3]*self.num_cones()

[docs]    def is_dcp(self):
"""An exponential constraint is DCP if each argument is affine.
"""
return all(arg.is_affine() for arg in self.args)

def canonicalize(self):
"""Canonicalizes by converting expressions to LinOps.
"""
arg_objs = []
arg_constr = []
for arg in self.args:
arg_objs.append(arg.canonical_form[0])
arg_constr + arg.canonical_form[1]
return 0, [ExpCone(*arg_objs)] + arg_constr

@pu.lazyprop
def __ECOS_format(self):
return ([], format_elemwise([self.x, self.z, self.y]))

@pu.lazyprop
def __SCS_format(self):
return ([], format_elemwise([self.x, self.y, self.z]))

@pu.lazyprop
def __CVXOPT_format(self):
constraints = []
for i, var in enumerate(self.vars_):
if var.type is not VARIABLE:
lone_var = lu.create_var(var.shape)
constraints.append(lu.create_eq(lone_var, var))
self.vars_[i] = lone_var
return (constraints, [])

def _solver_hook(self, vars_=None, scaling=None):
"""A function used by CVXOPT's nonlinear solver.

Based on f(x,y,z) = y * log(y) + x - y * log(z).

Parameters
----------
vars_: A cvxopt dense matrix with values for (x,y,z).
scaling: A scaling for the Hessian.

Returns
-------
_solver_hook() returns the constraint shape and a feasible point.
_solver_hook(x) returns the function value and gradient at x.
_solver_hook(x, z) returns the function value, gradient,
and (z scaled) Hessian at x.
"""
import cvxopt  # Not necessary unless using cvxopt solver.
entries = int(np.prod(self.shape))
if vars_ is None:
x_init = entries*[0.0]
y_init = entries*[0.5]
z_init = entries*[1.0]
return entries, cvxopt.matrix(x_init + y_init + z_init)
# Unpack vars_
x = vars_[0:entries]
y = vars_[entries:2*entries]
z = vars_[2*entries:]
# Out of domain.
# TODO what if y == 0.0?
if min(y) <= 0.0 or min(z) <= 0.0:
return None
# Evaluate the function.
f = cvxopt.matrix(0., (entries, 1))
for i in range(entries):
f[i] = x[i] - y[i]*math.log(z[i]) + y[i]*math.log(y[i])
# Compute the gradient.
Df = cvxopt.matrix(0., (entries, 3*entries))
for i in range(entries):
Df[i, i] = 1.0
Df[i, entries+i] = math.log(y[i]) - math.log(z[i]) + 1.0
Df[i, 2*entries+i] = -y[i]/z[i]

if scaling is None:
return f, Df
# Compute the Hessian.
big_H = cvxopt.spmatrix(0, [], [], size=(3*entries, 3*entries))
for i in range(entries):
H = cvxopt.matrix([
[0.0, 0.0, 0.0],
[0.0, 1.0/y[i], -1.0/z[i]],
[0.0, -1.0/z[i], y[i]/(z[i]**2)],
])
big_H[i:3*entries:entries, i:3*entries:entries] = scaling[i]*H
return f, Df, big_H