# Source code for pyomo.core.kernel.piecewise_library.transforms_nd

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  ___________________________________________________________________________

"""
This module contains transformations for representing a
multi-variate piecewise linear function using a
mixed-interger problem formulation. Reference::

Mixed-Integer Models for Non-separable Piecewise Linear \
Optimization: Unifying framework and Extensions (Vielma, \
Nemhauser 2008)
"""

import logging
import collections

from pyomo.core.kernel.block import block
from pyomo.core.kernel.set_types import IntegerSet
from pyomo.core.kernel.variable import (variable,
variable_dict,
variable_tuple)
from pyomo.core.kernel.constraint import (linear_constraint,
constraint_list,
constraint_tuple)
from pyomo.core.kernel.expression import (expression,
expression_tuple)
import pyomo.core.kernel.piecewise_library.util

logger = logging.getLogger('pyomo.core')

registered_transforms = {}

[docs]def piecewise_nd(tri,
values,
input=None,
output=None,
bound='eq',
repn='cc'):
"""
Models a multi-variate piecewise linear function.

This function takes a D-dimensional triangulation and a
list of function values associated with the points of
the triangulation and transforms this input data into a
block of variables and constraints that enforce a
piecewise linear relationship between an D-dimensional
vector of input variable and a single output
variable. In the general case, this transformation
requires the use of discrete decision variables.

Args:
tri (scipy.spatial.Delaunay): A triangulation over
the discretized variable domain. Can be
generated using a list of variables using the
utility function :func:util.generate_delaunay.
Required attributes:

- points: An (npoints, D) shaped array listing
the D-dimensional coordinates of the
discretization points.
- simplices: An (nsimplices, D+1) shaped array
of integers specifying the D+1 indices of
the points vector that define each simplex
of the triangulation.
values (numpy.array): An (npoints,) shaped array of
the values of the piecewise function at each of
coordinates in the triangulation points array.
input: A D-length list of variables or expressions
bound as the inputs of the piecewise function.
output: The variable constrained to be the output of
the piecewise linear function.
bound (str): The type of bound to impose on the
output expression. Can be one of:

- 'lb': y <= f(x)
- 'eq': y  = f(x)
- 'ub': y >= f(x)
repn (str): The type of piecewise representation to
use. Can be one of:

- 'cc': convex combination

Returns:
TransformedPiecewiseLinearFunctionND: a block \
containing any new variables, constraints, and \
other components used by the piecewise \
representation
"""
transform = None
try:
transform = registered_transforms[repn]
except KeyError:
raise ValueError(
"Keyword assignment repn='%s' is not valid. "
"Must be one of: %s"
% (repn,
str(sorted(registered_transforms.keys()))))
assert transform is not None

func = PiecewiseLinearFunctionND(tri,
values)

return transform(func,
input=input,
output=output,
bound=bound)

[docs]class PiecewiseLinearFunctionND(object):
"""A multi-variate piecewise linear function

Multi-varite piecewise linear functions are defined by a
triangulation over a finite domain and a list of
function values associated with the points of the
triangulation.  The function value between points in the
triangulation is implied through linear interpolation.

Args:
tri (scipy.spatial.Delaunay): A triangulation over
the discretized variable domain. Can be
generated using a list of variables using the
utility function :func:util.generate_delaunay.
Required attributes:

- points: An (npoints, D) shaped array listing
the D-dimensional coordinates of the
discretization points.
- simplices: An (nsimplices, D+1) shaped array
of integers specifying the D+1 indices of
the points vector that define each simplex
of the triangulation.
values (numpy.array): An (npoints,) shaped array of
the values of the piecewise function at each of
coordinates in the triangulation points array.
"""
__slots__ = ("_tri", "_values")

def __init__(self,
tri,
values,
validate=True,
**kwds):
assert pyomo.core.kernel.piecewise_library.util.numpy_available
assert pyomo.core.kernel.piecewise_library.util.scipy_available
assert isinstance(tri,
pyomo.core.kernel.piecewise_library.\
util.scipy.spatial.Delaunay)
assert isinstance(values,
pyomo.core.kernel.piecewise_library.\
util.numpy.ndarray)
npoints, ndim = tri.points.shape
nsimplices, _ = tri.simplices.shape
assert tri.simplices.shape == ndim + 1
assert nsimplices > 0
assert npoints > 0
assert ndim > 0
self._tri = tri
self._values = values

def __getstate__(self):
"""Required for older versions of the pickle
protocol since this class uses __slots__"""
return {key:getattr(self, key) for key in self.__slots__}

def __setstate__(self, state):
"""Required for older versions of the pickle
protocol since this class uses __slots__"""
for key in state:
setattr(self, key, state[key])

@property
def triangulation(self):
"""The triangulation over the domain of this function"""
return self._tri

@property
def values(self):
"""The set of values used to defined this function"""
return self._values

[docs]    def __call__(self, x):
"""
Evaluates the piecewise linear function using
interpolation. This method supports vectorized
function calls as the interpolation process can be
expensive for high dimensional data.

For the case when a single point is provided, the
argument x should be a (D,) shaped numpy array or
list, where D is the dimension of points in the
triangulation.

For the vectorized case, the argument x should be
a (n,D)-shaped numpy array.
"""
assert isinstance(x, collections.Sized)
if isinstance(x, pyomo.core.kernel.piecewise_library.\
util.numpy.ndarray):
if x.shape != self._tri.points.shape[1:]:
multi = True
assert x.shape[1:] == self._tri.points.shape, \
"%s != %s" % (x.shape, self._tri.points.shape)
else:
multi = False
else:
multi = False
_, ndim = self._tri.points.shape
i = self._tri.find_simplex(x)
if multi:
Tinv = self._tri.transform[i,:ndim]
r = self._tri.transform[i,ndim]
b = pyomo.core.kernel.piecewise_library.util.\
numpy.einsum('ijk,ik->ij', Tinv, x-r)
b = pyomo.core.kernel.piecewise_library.util.\
numpy.c_[b, 1 - b.sum(axis=1)]
s = self._tri.simplices[i]
return (b*self._values[s]).sum(axis=1)
else:
b = self._tri.transform[i,:ndim,:ndim].dot(
x - self._tri.transform[i,ndim,:])
s = self._tri.simplices[i]
val = b.dot(self._values[s[:ndim]])
val += (1-b.sum())*self._values[s[ndim]]
return val

[docs]class TransformedPiecewiseLinearFunctionND(block):
"""Base class for transformed multi-variate piecewise
linear functions

A transformed multi-variate piecewise linear functions
is a block of variables and constraints that enforce a
piecewise linear relationship between an vector input
variables and a single output variable.

Args:
f (:class:PiecewiseLinearFunctionND): The
multi-variate piecewise linear function to
transform.
input: The variable constrained to be the input of
the piecewise linear function.
output: The variable constrained to be the output of
the piecewise linear function.
bound (str): The type of bound to impose on the
output expression. Can be one of:

- 'lb': y <= f(x)
- 'eq': y  = f(x)
- 'ub': y >= f(x)
"""

def __init__(self,
f,
input=None,
output=None,
bound='eq'):
super(TransformedPiecewiseLinearFunctionND, self).__init__()
assert isinstance(f, PiecewiseLinearFunctionND)
if bound not in ('lb', 'ub', 'eq'):
raise ValueError("Invalid bound type %r. Must be "
"one of: ['lb','ub','eq']"
% (bound))
self._bound = bound
self._f = f
_,ndim = f._tri.points.shape
if input is None:
input = [None]*ndim
self._input = expression_tuple(
expression(input[i]) for i in range(ndim))
self._output = expression(output)

@property
def input(self):
"""The tuple of expressions that store the
inputs to the piecewise function. The returned
objects can be updated by assigning to their
:attr:expr attribute."""
return self._input

@property
def output(self):
"""The expression that stores the output of the
piecewise function. The returned object can be
updated by assigning to its :attr:expr
attribute."""
return self._output

@property
def bound(self):
"""The bound type assigned to the piecewise
relationship ('lb','ub','eq')."""
return self._bound

@property
def triangulation(self):
"""The triangulation over the domain of this function"""
return self._f.triangulation

@property
def values(self):
"""The set of values used to defined this function"""
return self._f.values

[docs]    def __call__(self, x):
"""
Evaluates the piecewise linear function using
interpolation. This method supports vectorized
function calls as the interpolation process can be
expensive for high dimensional data.

For the case when a single point is provided, the
argument x should be a (D,) shaped numpy array or
list, where D is the dimension of points in the
triangulation.

For the vectorized case, the argument x should be
a (n,D)-shaped numpy array.
"""
return self._f(x)

[docs]class piecewise_nd_cc(TransformedPiecewiseLinearFunctionND):
"""Discrete CC multi-variate piecewise representation

Expresses a multi-variate piecewise linear function
using the CC formulation.
"""

def __init__(self, *args, **kwds):
super(piecewise_nd_cc, self).__init__(*args, **kwds)

ndim = len(self.input)
nsimplices = len(self.triangulation.simplices)
npoints = len(self.triangulation.points)
pointsT = list(zip(*self.triangulation.points))

# create index objects
dimensions = range(ndim)
simplices = range(nsimplices)
vertices = range(npoints)

# create vars
self.v = variable_dict()
lmbda = self.v['lambda'] = variable_tuple(
variable(lb=0) for v in vertices)
y = self.v['y'] = variable_tuple(
variable(domain_type=IntegerSet, lb=0, ub=1) for s in simplices)
lmbda_tuple = tuple(lmbda)

# create constraints
self.c = constraint_list()

clist = []
for d in dimensions:
clist.append(linear_constraint(
variables=lmbda_tuple + (self.input[d],),
coefficients=tuple(pointsT[d]) + (-1,),
rhs=0))
self.c.append(constraint_tuple(clist))
del clist

self.c.append(linear_constraint(
variables=lmbda_tuple + (self.output,),
coefficients=tuple(self.values) + (-1,)))
if self.bound == 'ub':
self.c[-1].lb = 0
elif self.bound == 'lb':
self.c[-1].ub = 0
else:
assert self.bound == 'eq'
self.c[-1].rhs = 0

self.c.append(linear_constraint(
variables=lmbda_tuple,
coefficients=(1,)*len(lmbda_tuple),
rhs=1))

# generate a map from vertex index to simplex index,
# which avoids an n^2 lookup when generating the
# constraint
vertex_to_simplex = [[] for v in vertices]
for s, simplex in enumerate(self.triangulation.simplices):
for v in simplex:
vertex_to_simplex[v].append(s)

clist = []
for v in vertices:
variables = tuple(y[s] for s in vertex_to_simplex[v])
clist.append(linear_constraint(
variables=variables + (lmbda[v],),
coefficients=(1,)*len(variables) + (-1,),
lb=0))
self.c.append(constraint_tuple(clist))
del clist

self.c.append(linear_constraint(
variables=y,
coefficients=(1,)*len(y),
rhs=1))

registered_transforms['cc'] = piecewise_nd_cc