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/ProductPolynomialGradient.hxx"
22 #include "openturns/OSS.hxx"
23 #include "openturns/PersistentObjectFactory.hxx"
24
25 BEGIN_NAMESPACE_OPENTURNS
26
27
28
29 CLASSNAMEINIT(ProductPolynomialGradient)
30
31 static const Factory<ProductPolynomialGradient> Factory_ProductPolynomialGradient;
32
33
34 /* Default constructor */
ProductPolynomialGradient()35 ProductPolynomialGradient::ProductPolynomialGradient()
36 : GradientImplementation()
37 , polynomials_()
38 {
39 // Nothing to do
40 }
41
42
43 /* Constructor */
ProductPolynomialGradient(const PolynomialCollection & coll)44 ProductPolynomialGradient::ProductPolynomialGradient(const PolynomialCollection & coll)
45 : GradientImplementation()
46 , polynomials_(coll)
47 {
48 // Nothing to do
49 }
50
51
52 /* Virtual constructor */
clone() const53 ProductPolynomialGradient * ProductPolynomialGradient::clone() const
54 {
55 return new ProductPolynomialGradient(*this);
56 }
57
58
59 /* String converter */
__repr__() const60 String ProductPolynomialGradient::__repr__() const
61 {
62 return OSS() << "class=" << GetClassName();
63 }
64
65
66 /* Compute the gradient of a product of univariate polynomials */
gradient(const Point & inP) const67 Matrix ProductPolynomialGradient::gradient (const Point & inP) const
68 {
69 const UnsignedInteger inDimension = inP.getDimension();
70 if (inDimension != getInputDimension()) throw InvalidArgumentException(HERE) << "Error: trying to evaluate a ProductPolynomialFunction with an argument of invalid dimension";
71 Scalar productEvaluation = 1.0;
72 Point evaluations(inDimension);
73 Point derivatives(inDimension);
74 for (UnsignedInteger i = 0; i < inDimension; ++i)
75 {
76 const Scalar x = inP[i];
77 const Scalar y = polynomials_[i](x);
78 const Scalar dy = polynomials_[i].gradient(x);
79 evaluations[i] = y;
80 derivatives[i] = dy;
81 productEvaluation *= y;
82 }
83 Matrix grad(inDimension, 1);
84 // Usual case: productEvaluation <> 0
85 if (productEvaluation != 0.0)
86 {
87 for (UnsignedInteger i = 0; i < inDimension; ++i)
88 {
89 grad(i, 0) = derivatives[i] * (productEvaluation / evaluations[i]);
90 }
91 }
92 // Must compute the gradient in a more expensive way
93 else
94 {
95 for (UnsignedInteger i = 0; i < inDimension; ++i)
96 {
97 grad(i, 0) = derivatives[i];
98 for (UnsignedInteger j = 0; j < i; ++j) grad(i, 0) *= evaluations[j];
99 for (UnsignedInteger j = i + 1; j < inDimension; ++j) grad(i, 0) *= evaluations[j];
100 }
101 }
102 return grad;
103 }
104
105
106 /* Accessor for input point dimension */
getInputDimension() const107 UnsignedInteger ProductPolynomialGradient::getInputDimension() const
108 {
109 return polynomials_.getSize();
110 }
111
112 /* Accessor for output point dimension */
getOutputDimension() const113 UnsignedInteger ProductPolynomialGradient::getOutputDimension() const
114 {
115 return 1;
116 }
117
118
119 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const120 void ProductPolynomialGradient::save(Advocate & adv) const
121 {
122 GradientImplementation::save(adv);
123 adv.saveAttribute( "polynomials_", polynomials_ );
124 }
125
126 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)127 void ProductPolynomialGradient::load(Advocate & adv)
128 {
129 GradientImplementation::load(adv);
130 adv.loadAttribute( "polynomials_", polynomials_ );
131 }
132
133
134 END_NAMESPACE_OPENTURNS
135