# Source code for gemseo.core.mdofunctions.taylor_polynomials

# 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
#
# 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.
"""Functions computing first- and second-order Taylor polynomials from a function."""
from __future__ import annotations

from typing import Sequence

from numpy import ndarray

from gemseo.core.mdofunctions.mdo_function import ArrayType
from gemseo.core.mdofunctions.mdo_function import MDOFunction
from gemseo.core.mdofunctions.mdo_linear_function import MDOLinearFunction

[docs]def compute_linear_approximation(
function: MDOFunction,
x_vect: ArrayType,
name: str | None = None,
f_type: str | None = None,
input_names: Sequence[str] | None = None,
) -> MDOLinearFunction:
r"""Compute a first-order Taylor polynomial of a function.

:math:\newcommand{\xref}{\hat{x}}\newcommand{\dim}{d}
The first-order Taylor polynomial of a (possibly vector-valued) function
:math:f at a point :math:\xref is defined as

.. math::

\newcommand{\partialder}{\frac{\partial f}{\partial x_i}(\xref)}
f(x)
\approx
f(\xref) + \sum_{i = 1}^{\dim} \partialder \, (x_i - \xref_i).

Args:
function: The function to be linearized.
x_vect: The input vector at which to build the Taylor polynomial.
name: The name of the linear approximation function.
If None, create a name from the name of the function.
f_type: The type of the linear approximation function.
If None, the function will have no type.
input_names: The names of the inputs of the linear approximation function,
or a name base.
If None, use the names of the inputs of the function.

Returns:
The first-order Taylor polynomial of the function at the input vector.

Raises:
AttributeError: If the function does not have a Jacobian function.
"""
if not function.has_jac:
raise AttributeError("Function Jacobian unavailable for linear approximation.")

coefficients = function.jac(x_vect)
func_val = function.evaluate(x_vect)
if isinstance(func_val, ndarray):
# Make sure the function value is at most 1-dimensional
func_val = func_val.flatten()

return MDOLinearFunction(
coefficients,
f"{function.name}_linearized" if name is None else name,
f_type,
input_names if input_names else function.input_names,
func_val - coefficients @ x_vect,
)

function: MDOFunction,
x_vect: ArrayType,
hessian_approx: ArrayType,
input_names: Sequence[str] | None = None,
r"""Build a quadratic approximation of a function at a given point.

The function must be scalar-valued.

:math:\newcommand{\xref}{\hat{x}}\newcommand{\dim}{d}\newcommand{
\hessapprox}{\hat{H}}
For a given approximation :math:\hessapprox of the Hessian matrix of a
function :math:f at a point :math:\xref, the quadratic approximation of
:math:f is defined as

.. math::

\newcommand{\partialder}{\frac{\partial f}{\partial x_i}(\xref)}
f(x)
\approx
f(\xref)
+ \sum_{i = 1}^{\dim} \partialder \, (x_i - \xref_i)
+ \frac{1}{2} \sum_{i = 1}^{\dim} \sum_{j = 1}^{\dim}
\hessapprox_{ij} (x_i - \xref_i) (x_j - \xref_j).

Args:
function: The function to be approximated.
x_vect: The input vector at which to build the quadratic approximation.
hessian_approx: The approximation of the Hessian matrix
at this input vector.
input_names: The names of the inputs of the quadratic approximation function,
or a base name.
If None, use the ones of the current function.

Returns:
The second-order Taylor polynomial of the function at the given point.
"""
if (
not isinstance(hessian_approx, ndarray)
or hessian_approx.ndim != 2
or hessian_approx.shape != hessian_approx.shape
):
raise ValueError("Hessian approximation must be a square ndarray.")

if hessian_approx.shape != x_vect.size:
raise ValueError("Hessian approximation and vector must have same dimension.")

if not function.has_jac:
raise AttributeError("Jacobian unavailable.")