1 // -*- C++ -*-
2 /**
3 * @brief This is a nD polynomial build as a product of n 1D polynomial
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/ProductPolynomialHessian.hxx"
22 #include "openturns/OSS.hxx"
23 #include "openturns/PersistentObjectFactory.hxx"
24
25 BEGIN_NAMESPACE_OPENTURNS
26
27 CLASSNAMEINIT(ProductPolynomialHessian)
28
29 static const Factory<ProductPolynomialHessian> Factory_ProductPolynomialHessian;
30
31
32 /* Default constructor */
ProductPolynomialHessian()33 ProductPolynomialHessian::ProductPolynomialHessian()
34 : HessianImplementation()
35 , polynomials_()
36 {
37 // Nothing to do
38 }
39
40
41 /* Constructor */
ProductPolynomialHessian(const PolynomialCollection & coll)42 ProductPolynomialHessian::ProductPolynomialHessian(const PolynomialCollection & coll)
43 : HessianImplementation()
44 , polynomials_(coll)
45 {
46 // Nothing to do
47 }
48
49
50 /* Virtual constructor */
clone() const51 ProductPolynomialHessian * ProductPolynomialHessian::clone() const
52 {
53 return new ProductPolynomialHessian(*this);
54 }
55
56
57 /* String converter */
__repr__() const58 String ProductPolynomialHessian::__repr__() const
59 {
60 return OSS() << "class=" << GetClassName();
61 }
62
63
64 /* Compute the hessian of a product of univariate polynomials */
hessian(const Point & inP) const65 SymmetricTensor ProductPolynomialHessian::hessian (const Point & inP) const
66 {
67 const UnsignedInteger inDimension = inP.getDimension();
68 if (inDimension != getInputDimension()) throw InvalidArgumentException(HERE) << "Error: trying to evaluate a ProductPolynomialFunction with an argument of invalid dimension";
69 Scalar productEvaluation = 1.0;
70 Point evaluations(inDimension);
71 Point derivatives(inDimension);
72 Point secondDerivatives(inDimension);
73 for (UnsignedInteger i = 0; i < inDimension; ++i)
74 {
75 const Scalar x = inP[i];
76 const Scalar y = polynomials_[i](x);
77 const Scalar dy = polynomials_[i].gradient(x);
78 const Scalar d2y = polynomials_[i].hessian(x);
79 evaluations[i] = y;
80 derivatives[i] = dy;
81 secondDerivatives[i] = d2y;
82 productEvaluation *= y;
83 }
84 SymmetricTensor hess(inDimension, 1);
85 // Usual case: productEvaluation <> 0
86 if (productEvaluation != 0.0)
87 {
88 for (UnsignedInteger i = 0; i < inDimension; ++i)
89 {
90 const Scalar dyi = derivatives[i] * (productEvaluation / evaluations[i]);
91 for (UnsignedInteger j = 0; j < i; ++j)
92 {
93 hess(i, j, 0) = derivatives[j] * (dyi / evaluations[j]);
94 }
95 hess(i, i, 0) = secondDerivatives[i] * (productEvaluation / evaluations[i]);
96 }
97 }
98 // Must compute the hessian in a more expensive way
99 else
100 {
101 for (UnsignedInteger i = 0; i < inDimension; ++i)
102 {
103 for (UnsignedInteger j = 0; j < i; ++j)
104 {
105 hess(i, j, 0) = derivatives[i] * derivatives[j];
106 for (UnsignedInteger k = 0; k < j; ++k) hess(i, j, 0) *= evaluations[k];
107 for (UnsignedInteger k = j + 1; k < i; ++k) hess(i, j, 0) *= evaluations[k];
108 for (UnsignedInteger k = i + 1; k < inDimension; ++k) hess(i, j, 0) *= evaluations[k];
109 }
110 hess(i, i, 0) = secondDerivatives[i];
111 for (UnsignedInteger k = 0; k < i; ++k) hess(i, i, 0) *= evaluations[k];
112 for (UnsignedInteger k = i + 1; k < inDimension; ++k) hess(i, i, 0) *= evaluations[k];
113 }
114 }
115 return hess;
116 }
117
118 /* Accessor for input point dimension */
getInputDimension() const119 UnsignedInteger ProductPolynomialHessian::getInputDimension() const
120 {
121 return polynomials_.getSize();
122 }
123
124 /* Accessor for output point dimension */
getOutputDimension() const125 UnsignedInteger ProductPolynomialHessian::getOutputDimension() const
126 {
127 return 1;
128 }
129
130
131 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const132 void ProductPolynomialHessian::save(Advocate & adv) const
133 {
134 HessianImplementation::save(adv);
135 adv.saveAttribute( "polynomials_", polynomials_ );
136 }
137
138 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)139 void ProductPolynomialHessian::load(Advocate & adv)
140 {
141 HessianImplementation::load(adv);
142 adv.loadAttribute( "polynomials_", polynomials_ );
143 }
144
145
146 END_NAMESPACE_OPENTURNS
147