# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# This program 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
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# Contributors:
# INITIAL AUTHORS - API and implementation and/or documentation
# :author: Francois Gallard
# OTHER AUTHORS - MACROSCOPIC CHANGES
"""Constraints aggregation core functions."""
from __future__ import annotations
from math import log
from typing import TYPE_CHECKING
from numpy import argmax as np_argmax
from numpy import array
from numpy import atleast_2d
from numpy import exp as np_exp
from numpy import full
from numpy import heaviside
from numpy import max as np_max
from numpy import multiply
from numpy import ndarray
from numpy import sum as np_sum
from numpy import zeros
if TYPE_CHECKING:
from collections.abc import Sequence
[docs]
def compute_lower_bound_ks_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> float:
"""Transform a vector of constraint functions into a KS function.
The lower bound Kreisselmeier-Steinhauser function tends to the maximum operator
when the aggregation parameter tends to infinity.
See :cite:`kennedy2015improved` and :cite:`kreisselmeier1983application`.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The aggregation parameter.
scale: The scaling factor for multiplying the constraints.
Returns:
The KS function value.
"""
alpha = len(orig_val)
return compute_upper_bound_ks_agg(orig_val, indices, rho, scale) - log(alpha) / rho
[docs]
def compute_upper_bound_ks_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> float:
"""Transform a vector of constraint functions into a KS function.
The upper bound Kreisselmeier-Steinhauser function tends to the maximum operator
when the aggregation parameter tends to infinity.
See :cite:`kennedy2015improved` and :cite:`kreisselmeier1983application`.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The aggregation parameter.
scale: The scaling factor for multiplying the constraints.
Returns:
The KS function value.
"""
if indices is not None:
orig_val = orig_val[indices]
orig_val *= scale
m = max(orig_val)
return m + (1.0 / rho) * log(sum(np_exp(rho * (orig_val + 1.0 - m)))) - 1.0
[docs]
def compute_total_ks_agg_jac(
orig_val: ndarray,
orig_jac: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of KS function with respect to constraint function inputs.
See :cite:`kennedy2015improved` and :cite:`kreisselmeier1983application`.
Args:
orig_val: The original constraint values.
orig_jac: The original constraint jacobian.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The multiplicative parameter in the exponential.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of KS function with respect to constraint function inputs.
"""
if indices is not None:
orig_jac = orig_jac[indices, :]
orig_val = orig_val[indices]
orig_jac *= scale
orig_val *= scale
m = max(orig_val)
div = np_sum(np_exp(rho * (orig_val + 1.0 - m)))
weights = np_exp(rho * (orig_val + 1.0 - m)).T / div
return np_sum(multiply(atleast_2d(weights).T, orig_jac), axis=0)
[docs]
def compute_partial_ks_agg_jac(
orig_val: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of KS function with respect to constraint functions.
See :cite:`kennedy2015improved` and :cite:`kreisselmeier1983application`.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The multiplicative parameter in the exponential.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of KS function with respect to constraint functions.
"""
full_size = orig_val.size
if indices is not None:
orig_val = orig_val[indices]
orig_val *= scale
m = max(orig_val)
div = np_sum(np_exp(rho * (orig_val + 1.0 - m)))
weights = np_exp(rho * (orig_val + 1.0 - m)).T / div
der = atleast_2d(multiply(weights, scale))
return __filter_jac(der, full_size, indices)
[docs]
def compute_iks_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> float:
"""Transform a vector of constraint functions into an induces exponential function.
The induces exponential function (IKS) tends to the maximum operator
when the aggregation parameter tends to infinity.
See :cite:`kennedy2015improved`.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The multiplicative parameter in the exponential.
scale: The scaling factor for multiplying the constraints.
Returns:
The IKS function value.
"""
if indices is not None:
orig_val = orig_val[indices]
orig_val *= scale
m = max(orig_val)
iks = sum(orig_val * np_exp(rho * (orig_val + 1.0 - m)))
iks /= sum(np_exp(rho * (orig_val + 1.0 - m)))
return iks
[docs]
def compute_total_iks_agg_jac(
orig_val: ndarray,
orig_jac: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of IKS function with respect to constraints inputs.
See :cite:`kennedy2015improved`.
Args:
orig_val: The original constraint values.
orig_jac: The original constraint jacobian.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The multiplicative parameter in the exponential.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of IKS function with respect to constraints inputs.
"""
if indices is not None:
orig_jac = orig_jac[indices, :]
orig_val = orig_val[indices]
orig_jac *= scale
orig_val *= scale
m = max(orig_val)
iks_num = sum(orig_val * np_exp(rho * (orig_val + 1.0 - m)))
iks_den = sum(np_exp(rho * (orig_val + 1.0 - m)))
iks_num_der = np_sum(
multiply(atleast_2d(np_exp(rho * (orig_val + 1.0 - m))).T, orig_jac), axis=0
)
iks_num_der += np_sum(
multiply(
atleast_2d(np_exp(rho * (orig_val + 1.0 - m)) * orig_val).T,
rho * orig_jac,
),
axis=0,
)
iks_den_der = np_sum(
multiply(atleast_2d(np_exp(rho * (orig_val + 1.0 - m))).T, rho * orig_jac),
axis=0,
)
return (-iks_den_der / iks_den**2) * iks_num + iks_num_der / iks_den
[docs]
def compute_partial_iks_agg_jac(
orig_val: ndarray,
indices: Sequence[int] | None = None,
rho: float = 1e2,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of IKS function with respect to constraints functions.
See :cite:`kennedy2015improved`.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
rho: The multiplicative parameter in the exponential.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of IKS function with respect to constraints functions.
"""
full_size = orig_val.size
if indices is not None:
orig_val = orig_val[indices]
orig_val *= scale
m = max(orig_val)
iks_num = sum(orig_val * np_exp(rho * (orig_val + 1.0 - m)))
iks_den = sum(np_exp(rho * (orig_val + 1.0 - m)))
iks_num_der = atleast_2d(np_exp(rho * (orig_val + 1.0 - m)))
iks_num_der += rho * atleast_2d(np_exp(rho * (orig_val + 1.0 - m)) * orig_val)
iks_den_der = rho * atleast_2d(np_exp(rho * (orig_val + 1.0 - m)))
iks_d = atleast_2d((-iks_den_der / iks_den**2) * iks_num + iks_num_der / iks_den)
iks_d = multiply(iks_d, scale)
return __filter_jac(iks_d, full_size, indices)
def __filter_jac(
orig_jac: ndarray,
full_size: int,
indices: Sequence[int],
) -> ndarray:
"""Filters the Jacobian according to the indices.
Args:
orig_jac: The original Jacobian.
full_size: The size of the full flatten Jacobian array.
indices: The array of indices to filter.
Returns:
The filtered Jacobian.
"""
if indices is not None:
jac = zeros((1, full_size))
jac[:, indices] = orig_jac
return jac
return orig_jac
[docs]
def compute_sum_square_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Transform a vector of constraint functions into a sum squared function.
The sum squared function is the sum of the squares of the input vector components.
Args:
orig_val: The input vector.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, scale all the constraint values.
scale: The scaling factor for multiplying the constraints.
Returns:
The sum squared function value.
"""
if indices is not None:
orig_val = orig_val[indices]
return np_sum(scale * orig_val**2)
[docs]
def compute_total_sum_square_agg_jac(
orig_val: ndarray,
orig_jac: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of squared sum function with respect to constraints inputs.
Args:
orig_val: The original constraint values.
orig_jac: The original constraint jacobian.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of squared sum function with respect to constraints inputs.
"""
orig_jac = atleast_2d(orig_jac)
if indices is not None:
orig_jac = orig_jac[indices, :]
orig_val = orig_val[indices]
return np_sum((2 * scale * orig_val).flatten() * orig_jac.T, axis=1)
[docs]
def compute_partial_sum_square_agg_jac(
orig_val: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of squared sum function with respect to constraints.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of squared sum function with respect to constraints.
"""
if indices is not None:
jac = zeros((1, orig_val.size))
jac[:, indices] = 2.0 * scale * orig_val[indices]
else:
jac = full((1, orig_val.size), 2.0) * atleast_2d(scale * orig_val)
return jac
[docs]
def compute_max_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> float:
"""Transform a vector of constraints into a max of all values.
The maximum function is not differentiable for all input values.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The maximum value.
"""
if indices is not None:
orig_val = orig_val[indices]
orig_val *= scale
return array([np_max(orig_val)])
[docs]
def compute_max_agg_jac(
orig_val: ndarray,
orig_jac: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of max function with respect to constraints inputs.
Args:
orig_val: The original constraint values.
orig_jac: The original constraint jacobian.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of max function with respect to constraints inputs.
"""
if indices is not None:
orig_jac = orig_jac[indices, :]
orig_val = orig_val[indices]
orig_jac *= scale
orig_val *= scale
i_max = np_argmax(orig_val)
return atleast_2d(orig_jac)[i_max, :]
[docs]
def compute_sum_positive_square_agg(
orig_val: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Transform a vector of constraint functions into a positive sum squared function.
The positive sum squared function is the sum of the squares of the input vector
components that are positive.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, scale all the constraint values.
scale: The scaling factor for multiplying the constraints.
Returns:
The positive sum squared function value.
"""
if indices is not None:
orig_val = orig_val[indices]
return np_sum(scale * (orig_val**2) * heaviside(orig_val.real, 0))
[docs]
def compute_total_sum_square_positive_agg_jac(
orig_val: ndarray,
orig_jac: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of positive sum squared function w.r.t. constraints inputs.
Args:
orig_val: The original constraint values.
orig_jac: The original constraint jacobian.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of positive sum squared function w.r.t. constraints inputs.
"""
orig_jac = atleast_2d(orig_jac)
if indices is not None:
orig_jac = orig_jac[indices, :]
orig_val = orig_val[indices]
return np_sum(
(2 * scale * orig_val * heaviside(orig_val.real, 0)).flatten() * orig_jac.T,
axis=1,
)
[docs]
def compute_partial_sum_positive_square_agg_jac(
orig_val: ndarray,
indices: Sequence[int] | None = None,
scale: float | ndarray = 1.0,
) -> ndarray:
"""Compute the Jacobian of positive sum squared function w.r.t. constraints.
Args:
orig_val: The original constraint values.
indices: The indices to generate a subset of the outputs to aggregate.
If ``None``, aggregate all the outputs.
scale: The scaling factor for multiplying the constraints.
Returns:
The Jacobian of positive sum squared function w.r.t. constraints.
"""
if indices is not None:
jac = zeros((1, orig_val.size))
jac[:, indices] = (
2.0 * scale * orig_val[indices] * heaviside(orig_val[indices], 0.0)
)
else:
jac = atleast_2d(2 * scale * orig_val * heaviside(orig_val.real, 0.0))
return jac
