Source code for gemseo.utils.derivatives.complex_step

# -*- coding: utf-8 -*-
# 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
"""Gradient approximation by complex step."""
from __future__ import division, unicode_literals

import logging
from typing import Any, Callable, List, Optional, Sequence, Tuple, Union

from numpy import complex128, finfo, ndarray, where, zeros
from numpy.linalg import norm

from gemseo.algos.design_space import DesignSpace
from gemseo.core.parallel_execution import ParallelExecution
from gemseo.utils.derivatives.gradient_approximator import GradientApproximator
from gemseo.utils.py23_compat import xrange

EPSILON = finfo(float).eps
LOGGER = logging.getLogger(__name__)


[docs]class ComplexStep(GradientApproximator): r"""Complex step approximator, performing a second-order gradient calculation. Enable a much lower step than real finite differences, typically 1e-30, since there is no cancellation error due to a difference calculation. .. math:: \frac{df(x)}{dx} \approx Im\left( \frac{f(x+j*\\delta x)}{\\delta x} \right) See Martins, Joaquim RRA, Peter Sturdza, and Juan J. Alonso. "The complex-step derivative approximation." ACM Transactions on Mathematical Software (TOMS) 29.3 (2003): 245-262. """ ALIAS = "complex_step" def __init__( self, f_pointer, # type: Callable[[ndarray],ndarray] step=1e-20, # type: complex parallel=False, # type: bool design_space=None, # type: Optional[DesignSpace] normalize=True, # type: bool **parallel_args # type: Union[int,bool,float] ): # type: (...) -> None super(ComplexStep, self).__init__( f_pointer, step=step, parallel=parallel, design_space=design_space, normalize=True, **parallel_args ) @GradientApproximator.step.setter def step(self, value): if value.imag != 0: self._step = value.imag else: self._step = value
[docs] def f_gradient( self, x_vect, # type: ndarray step=None, # type: Optional[complex] x_indices=None, # type: Optional[Sequence[int]] **kwargs # type: Any ): # type: (...) -> ndarray if norm(x_vect.imag) != 0.0: raise ValueError( "Impossible to check the gradient at a complex " "point using the complex step method." ) return super(ComplexStep, self).f_gradient( x_vect, step=step, x_indices=x_indices, **kwargs )
def _compute_parallel_grad( self, input_values, # type: ndarray n_perturbations, # type: int input_perturbations, # type: ndarray step, # type: float **kwargs # type: Any ): # type: (...) -> ndarray def func_noargs( f_input_values, # type: ndarray ): # type:(...) -> ndarray """Call the function without explicitly passed arguments. Args: f_input_values: The input value. Return: The value of the function output. """ return self.f_pointer(f_input_values, **kwargs) functions = [func_noargs] * n_perturbations parallel_execution = ParallelExecution(functions, **self._par_args) perturbated_inputs = [ input_values + input_perturbations[:, perturbation_index] for perturbation_index in xrange(n_perturbations) ] perturbated_outputs = parallel_execution.execute(perturbated_inputs) gradient = [] for perturbation_index in xrange(n_perturbations): gradient.append( perturbated_outputs[perturbation_index].imag / input_perturbations[perturbation_index, perturbation_index].imag ) return gradient def _compute_grad( self, input_values, # type: ndarray n_perturbations, # type: int input_perturbations, # type: ndarray step, # type: float **kwargs # type: Any ): # type: (...) -> ndarray gradient = [] for perturbation_index in xrange(n_perturbations): perturbated_input = ( input_values + input_perturbations[:, perturbation_index] ) perturbated_output = self.f_pointer(perturbated_input, **kwargs) gradient.append( perturbated_output.imag / input_perturbations[perturbation_index, perturbation_index].imag ) return gradient def _generate_perturbations( self, input_values, # type: ndarray input_indices, # type: List[int] step, # type: float ): # type: (...) -> Tuple[ndarray,Union[float,ndarray]] input_dimension = len(input_values) n_indices = len(input_indices) input_perturbations = zeros((input_dimension, n_indices), dtype=complex128) x_nnz = where(input_values == 0.0, 1.0, input_values)[input_indices] input_perturbations[input_indices, range(n_indices)] = 1j * x_nnz * step return input_perturbations, step