1 // -*- C++ -*-
2 /**
3 * @brief This is a nD function build as a product of n 1D function
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include "openturns/ProductUniVariateFunctionHessian.hxx"
22 #include "openturns/OSS.hxx"
23 #include "openturns/PersistentObjectFactory.hxx"
24
25 BEGIN_NAMESPACE_OPENTURNS
26
27 CLASSNAMEINIT(ProductUniVariateFunctionHessian)
28
29 static const Factory<ProductUniVariateFunctionHessian> Factory_ProductUniVariateFunctionHessian;
30
31
32 /* Default constructor */
ProductUniVariateFunctionHessian()33 ProductUniVariateFunctionHessian::ProductUniVariateFunctionHessian()
34 : HessianImplementation()
35 {
36 // Nothing to do
37 }
38
39
40 /* Constructor */
ProductUniVariateFunctionHessian(const ProductUniVariateFunctionEvaluation & evaluation)41 ProductUniVariateFunctionHessian::ProductUniVariateFunctionHessian(const ProductUniVariateFunctionEvaluation & evaluation)
42 : HessianImplementation()
43 , p_evaluation_(evaluation.clone())
44 {
45 // Nothing to do
46 }
47
48
49 /* Constructor */
ProductUniVariateFunctionHessian(const Pointer<ProductUniVariateFunctionEvaluation> & p_evaluation)50 ProductUniVariateFunctionHessian::ProductUniVariateFunctionHessian(const Pointer<ProductUniVariateFunctionEvaluation> & p_evaluation)
51 : HessianImplementation()
52 , p_evaluation_(p_evaluation)
53 {
54 // Nothing to do
55 }
56
57
58 /* Virtual constructor */
clone() const59 ProductUniVariateFunctionHessian * ProductUniVariateFunctionHessian::clone() const
60 {
61 return new ProductUniVariateFunctionHessian(*this);
62 }
63
64
65 /* String converter */
__repr__() const66 String ProductUniVariateFunctionHessian::__repr__() const
67 {
68 return OSS() << "class=" << GetClassName();
69 }
70
71
72 /* Compute the hessian of a product of univariate polynomials */
hessian(const Point & inP) const73 SymmetricTensor ProductUniVariateFunctionHessian::hessian (const Point & inP) const
74 {
75 const UnsignedInteger inDimension = inP.getDimension();
76 if (inDimension != getInputDimension()) throw InvalidArgumentException(HERE) << "Error: trying to evaluate a ProductPolynomialFunction with an argument of invalid dimension";
77 Scalar productEvaluation = 1.0;
78 Point evaluations(inDimension);
79 Point derivatives(inDimension);
80 Point secondDerivatives(inDimension);
81 for (UnsignedInteger i = 0; i < inDimension; ++i)
82 {
83 const Scalar x = inP[i];
84 const Scalar y = p_evaluation_->functions_[i](x);
85 const Scalar dy = p_evaluation_->functions_[i].gradient(x);
86 const Scalar d2y = p_evaluation_->functions_[i].hessian(x);
87 evaluations[i] = y;
88 derivatives[i] = dy;
89 secondDerivatives[i] = d2y;
90 productEvaluation *= y;
91 }
92 SymmetricTensor hess(inDimension, 1);
93 // Usual case: productEvaluation <> 0
94 if (productEvaluation != 0.0)
95 {
96 for (UnsignedInteger i = 0; i < inDimension; ++i)
97 {
98 const Scalar dyi = derivatives[i] * (productEvaluation / evaluations[i]);
99 for (UnsignedInteger j = 0; j < i; ++j)
100 {
101 hess(i, j, 0) = derivatives[j] * (dyi / evaluations[j]);
102 }
103 hess(i, i, 0) = secondDerivatives[i] * (productEvaluation / evaluations[i]);
104 }
105 }
106 // Must compute the hessian in a more expensive way
107 else
108 {
109 for (UnsignedInteger i = 0; i < inDimension; ++i)
110 {
111 for (UnsignedInteger j = 0; j < i; ++j)
112 {
113 hess(i, j, 0) = derivatives[i] * derivatives[j];
114 for (UnsignedInteger k = 0; k < j; ++k) hess(i, j, 0) *= evaluations[k];
115 for (UnsignedInteger k = j + 1; k < i; ++k) hess(i, j, 0) *= evaluations[k];
116 for (UnsignedInteger k = i + 1; k < inDimension; ++k) hess(i, j, 0) *= evaluations[k];
117 }
118 hess(i, i, 0) = secondDerivatives[i];
119 for (UnsignedInteger k = 0; k < i; ++k) hess(i, i, 0) *= evaluations[k];
120 for (UnsignedInteger k = i + 1; k < inDimension; ++k) hess(i, i, 0) *= evaluations[k];
121 }
122 }
123 return hess;
124 }
125
126 /* Accessor for input point dimension */
getInputDimension() const127 UnsignedInteger ProductUniVariateFunctionHessian::getInputDimension() const
128 {
129 return p_evaluation_->functions_.getSize();
130 }
131
132 /* Accessor for output point dimension */
getOutputDimension() const133 UnsignedInteger ProductUniVariateFunctionHessian::getOutputDimension() const
134 {
135 return 1;
136 }
137
138
139 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const140 void ProductUniVariateFunctionHessian::save(Advocate & adv) const
141 {
142 HessianImplementation::save(adv);
143 adv.saveAttribute("evaluation_", *p_evaluation_);
144 }
145
146 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)147 void ProductUniVariateFunctionHessian::load(Advocate & adv)
148 {
149 HessianImplementation::load(adv);
150 TypedInterfaceObject<ProductUniVariateFunctionEvaluation> evaluation;
151 adv.loadAttribute("evaluation_", evaluation);
152 p_evaluation_ = evaluation.getImplementation();
153 }
154
155
156 END_NAMESPACE_OPENTURNS
157