1 #include "MUQ/Approximation/Polynomials/BasisExpansion.h"
2 #include "MUQ/Approximation/Polynomials/MonotoneExpansion.h"
3
4 #include "MUQ/Approximation/Polynomials/Monomial.h"
5 #include "MUQ/Utilities/MultiIndices/MultiIndexFactory.h"
6
7 #include "gtest/gtest.h"
8
9 using namespace muq::Approximation;
10 using namespace muq::Utilities;
11
12
13 class Approximation_MonotoneExpansion1d : public::testing::Test {
14 public:
Approximation_MonotoneExpansion1d()15 Approximation_MonotoneExpansion1d() {
16
17 auto monomial = std::make_shared<Monomial>();
18 auto bases = std::vector<std::shared_ptr<IndexedScalarBasis>>(1, monomial);
19
20 std::shared_ptr<MultiIndexSet> multis = MultiIndexFactory::CreateTotalOrder(1, 2);
21
22 monoPart = std::make_shared<BasisExpansion>(bases, multis);
23 expansion1 = std::make_shared<MonotoneExpansion>(monoPart);
24 expansion2 = std::make_shared<MonotoneExpansion>(monoPart, true);
25
26 coeffs.resize(4);
27 coeffs << 0.0, -0.01, 0.1, -0.5; // intercept, slope
28
29 expansion1->SetCoeffs(coeffs);
30 expansion2->SetCoeffs(coeffs);
31 }
32
~Approximation_MonotoneExpansion1d()33 virtual ~Approximation_MonotoneExpansion1d() {};
34
35 std::shared_ptr<BasisExpansion> monoPart;
36 std::shared_ptr<MonotoneExpansion> expansion1, expansion2;
37 Eigen::VectorXd coeffs;
38
39 };
40
TEST_F(Approximation_MonotoneExpansion1d,Evaluate)41 TEST_F(Approximation_MonotoneExpansion1d, Evaluate){
42
43 Eigen::VectorXd evalPt(1);
44
45 int numSteps = 10;
46 double ub = 2.0;
47 double lb = -2.0;
48 double dx = (ub-lb)/numSteps;
49
50 evalPt << lb;
51
52 double oldOutput = expansion1->Evaluate(evalPt).at(0)(0);
53
54 for(int i=1; i<numSteps; ++i){
55 evalPt(0) = lb + dx*double(i);
56 double newOutput = expansion1->Evaluate(evalPt).at(0)(0);
57 EXPECT_GT(newOutput, oldOutput);
58
59 double monoEval = monoPart->Evaluate(evalPt).at(0)(0);
60 double trueMono = coeffs(1) + coeffs(2)*evalPt(0) + coeffs(3)*evalPt(0)*evalPt(0);
61 EXPECT_DOUBLE_EQ(trueMono, monoEval);
62
63 double truth = coeffs(1)*coeffs(1)*evalPt(0)
64 + coeffs(1)*coeffs(2)*std::pow(evalPt(0),2.0)
65 + (1.0/3.0)*(2.0*coeffs(1)*coeffs(3)+coeffs(2)*coeffs(2))*std::pow(evalPt(0),3.0)
66 + 0.5*coeffs(2)*coeffs(3)*std::pow(evalPt(0),4.0)
67 + 0.2*coeffs(3)*coeffs(3)*std::pow(evalPt(0),5.0);
68 EXPECT_NEAR(truth, newOutput, 1e-4);
69
70
71 Eigen::MatrixXd jac = expansion1->Jacobian(0,0,evalPt);
72 EXPECT_GT(jac(0,0), 0.0);
73
74 oldOutput = newOutput;
75 }
76 }
77
78
TEST_F(Approximation_MonotoneExpansion1d,JacobianFD)79 TEST_F(Approximation_MonotoneExpansion1d, JacobianFD){
80
81 Eigen::VectorXd evalPt(1);
82
83 int numSteps = 10;
84 double ub = 2.0;
85 double lb = -2.0;
86 double dx = (ub-lb)/numSteps;
87
88 evalPt << lb;
89
90 const double eps = 1e-3;
91
92 for(int i=1; i<numSteps; ++i){
93
94 evalPt(0) = lb + dx*double(i);
95 double f1 = expansion1->Evaluate(evalPt).at(0)(0);
96
97 Eigen::MatrixXd jac = expansion1->Jacobian(0,0,evalPt);
98
99 evalPt(0) += eps;
100 double f2 = expansion1->Evaluate(evalPt).at(0)(0);
101
102 EXPECT_NEAR((f2-f1)/eps, jac(0,0), 3e-3);
103 }
104 }
105
TEST_F(Approximation_MonotoneExpansion1d,CoeffJacobian)106 TEST_F(Approximation_MonotoneExpansion1d, CoeffJacobian){
107
108 Eigen::VectorXd evalPt(1);
109 evalPt << 0.25;
110
111 const double eps = 1e-3;
112
113 double f1 = expansion1->Evaluate(evalPt).at(0)(0);
114 Eigen::MatrixXd jac = expansion2->Jacobian(1, 0, evalPt, coeffs);
115
116 for(int i=0; i<coeffs.size(); ++i){
117 Eigen::VectorXd newCoeffs = coeffs;
118 newCoeffs(i) += eps;
119
120 double f2 = expansion2->Evaluate(evalPt, newCoeffs).at(0)(0);
121
122 EXPECT_NEAR((f2-f1)/eps, jac(0,i), 1e-3);
123 }
124 }
125
TEST_F(Approximation_MonotoneExpansion1d,Determinant)126 TEST_F(Approximation_MonotoneExpansion1d, Determinant){
127
128 Eigen::VectorXd evalPt(1);
129 evalPt << 0.25;
130
131 double logDet1 = expansion1->LogDeterminant(evalPt);
132 Eigen::VectorXd grad = expansion2->GradLogDeterminant(evalPt,coeffs);
133
134 const double eps = 1e-7;
135
136 for(int i=0; i<coeffs.size(); ++i){
137 Eigen::VectorXd newCoeffs = coeffs;
138 newCoeffs(i) += eps;
139 double logDet2 = expansion2->LogDeterminant(evalPt,newCoeffs);
140
141 EXPECT_NEAR((logDet2-logDet1)/eps, grad(i), 1e-3);
142 }
143 }
144
145
TEST_F(Approximation_MonotoneExpansion1d,EvaluateWithCoeffs)146 TEST_F(Approximation_MonotoneExpansion1d, EvaluateWithCoeffs){
147
148 Eigen::VectorXd evalPt(1);
149
150 Eigen::VectorXd newCoeffs(4);
151 newCoeffs << 0.1, 0.5, 2.0, -0.1;
152 coeffs = newCoeffs;
153
154 int numSteps = 10;
155 double ub = 2.0;
156 double lb = -2.0;
157 double dx = (ub-lb)/numSteps;
158
159 evalPt << lb;
160 double oldOutput = expansion2->Evaluate(evalPt, newCoeffs).at(0)(0);
161
162 for(int i=1; i<numSteps; ++i){
163 evalPt(0) = lb + dx*double(i);
164 double newOutput = expansion2->Evaluate(evalPt, newCoeffs).at(0)(0);
165 EXPECT_GT(newOutput, oldOutput);
166
167 double truth = coeffs(0) + coeffs(1)*coeffs(1)*evalPt(0)
168 + coeffs(1)*coeffs(2)*std::pow(evalPt(0),2.0)
169 + (1.0/3.0)*(2.0*coeffs(1)*coeffs(3)+coeffs(2)*coeffs(2))*std::pow(evalPt(0),3.0)
170 + 0.5*coeffs(2)*coeffs(3)*std::pow(evalPt(0),4.0)
171 + 0.2*coeffs(3)*coeffs(3)*std::pow(evalPt(0),5.0);
172 EXPECT_NEAR(truth, newOutput, 5e-2);
173
174 oldOutput = newOutput;
175 }
176 }
177
178 class Approximation_MonotoneExpansion2d : public::testing::Test {
179 public:
Approximation_MonotoneExpansion2d()180 Approximation_MonotoneExpansion2d() {
181
182 auto monomial = std::make_shared<Monomial>();
183
184 // Build the general pieces
185 auto bases = std::vector<std::shared_ptr<IndexedScalarBasis>>(2, monomial);
186 std::shared_ptr<MultiIndexSet> constantMulti = MultiIndexFactory::CreateTotalOrder(2, 0);
187 std::shared_ptr<MultiIndexSet> multis = MultiIndexFactory::CreateTotalOrder(2, 2, 0, std::make_shared<DimensionLimiter>(0,1));
188 Eigen::MatrixXd coeffs(1,3);
189 coeffs << 0.0, 0.5, -0.1;
190 generalParts.push_back(std::make_shared<BasisExpansion>(bases, constantMulti));
191 generalParts.push_back(std::make_shared<BasisExpansion>(bases, multis, coeffs));
192
193 // Build the monotone pieces
194 coeffs << 0.1, 1.0, 0.5;
195 monoParts.push_back(std::make_shared<BasisExpansion>(bases, multis, coeffs));
196
197 multis = MultiIndexFactory::CreateTotalOrder(2, 2);
198 coeffs = Eigen::MatrixXd::Random(1,multis->Size());
199 monoParts.push_back(std::make_shared<BasisExpansion>(bases, multis, coeffs));
200
201 // Construct the expansion
202 expansion1 = std::make_shared<MonotoneExpansion>(generalParts, monoParts);
203 expansion2 = std::make_shared<MonotoneExpansion>(generalParts, monoParts, true);
204 };
205
~Approximation_MonotoneExpansion2d()206 virtual ~Approximation_MonotoneExpansion2d() {};
207
208 std::vector<std::shared_ptr<BasisExpansion>> generalParts;
209 std::vector<std::shared_ptr<BasisExpansion>> monoParts;
210
211 std::shared_ptr<MonotoneExpansion> expansion1, expansion2;
212 };
213
TEST_F(Approximation_MonotoneExpansion2d,Evaluate)214 TEST_F(Approximation_MonotoneExpansion2d, Evaluate){
215
216 Eigen::VectorXd evalPt = Eigen::VectorXd::Zero(2);
217
218 int numSteps = 10;
219 double ub = 2.0;
220 double lb = -2.0;
221 double dx = (ub-lb)/numSteps;
222
223 evalPt(1) = lb;
224 Eigen::VectorXd oldOutput = expansion1->Evaluate(evalPt).at(0);
225
226 for(int i=1; i<numSteps; ++i){
227 evalPt(1) = lb + dx*double(i);
228 Eigen::VectorXd newOutput = expansion1->Evaluate(evalPt).at(0);
229 EXPECT_GT(newOutput(1), oldOutput(1));
230
231 oldOutput = newOutput;
232 }
233
234 }
235
TEST_F(Approximation_MonotoneExpansion2d,Inverse)236 TEST_F(Approximation_MonotoneExpansion2d, Inverse){
237
238 Eigen::VectorXd tgtPt = Eigen::Vector2d{0.1, 0.5};
239 Eigen::VectorXd refPt = expansion1->EvaluateForward(tgtPt);
240 Eigen::VectorXd invPt = expansion1->EvaluateInverse(refPt);
241
242 EXPECT_NEAR(tgtPt(0), invPt(0), 1e-8);
243 EXPECT_NEAR(tgtPt(1), invPt(1), 1e-8);
244
245 tgtPt = Eigen::Vector2d{-9.2, 5.5};
246 refPt = expansion1->EvaluateForward(tgtPt);
247 invPt = expansion1->EvaluateInverse(refPt);
248
249 EXPECT_NEAR(tgtPt(0), invPt(0), 1e-8);
250 EXPECT_NEAR(tgtPt(1), invPt(1), 1e-8);
251 }
252
253
254
TEST_F(Approximation_MonotoneExpansion2d,Jacobian)255 TEST_F(Approximation_MonotoneExpansion2d, Jacobian){
256
257 Eigen::VectorXd evalPt(2);
258 evalPt << 0.1, 0.2;
259
260 Eigen::VectorXd f1 = expansion1->Evaluate(evalPt).at(0);
261 Eigen::MatrixXd jac = expansion1->Jacobian(0,0, evalPt);
262
263 EXPECT_GT(jac(0,0),0.0);
264 EXPECT_GT(jac(1,1),0.0);
265
266 const double eps = 1e-6;
267 Eigen::VectorXd newPt = evalPt;
268 newPt(0) += eps;
269 Eigen::VectorXd f2 = expansion1->Evaluate(newPt).at(0);
270
271 for(int i=0; i<f2.size(); ++i)
272 EXPECT_NEAR((f2(i)-f1(i))/eps, jac(i,0), 1e-3);
273
274 newPt = evalPt;
275 newPt(1) += eps;
276 f2 = expansion1->Evaluate(newPt).at(0);
277
278 for(int i=0; i<f2.size(); ++i)
279 EXPECT_NEAR((f2(i)-f1(i))/eps, jac(i,1), 1e-3);
280 }
281
TEST_F(Approximation_MonotoneExpansion2d,CoeffsJacobian)282 TEST_F(Approximation_MonotoneExpansion2d, CoeffsJacobian){
283
284 Eigen::VectorXd evalPt(2);
285 evalPt << 0.1, 0.2;
286
287 Eigen::VectorXd coeffs = expansion2->GetCoeffs();
288
289 const double eps = 1e-3;
290
291 Eigen::VectorXd f1 = expansion2->Evaluate(evalPt, coeffs).at(0);
292 Eigen::MatrixXd jac = expansion2->Jacobian(1,0,evalPt,coeffs);
293
294 for(int i=0; i<coeffs.size(); ++i){
295 Eigen::VectorXd newCoeffs = coeffs;
296 newCoeffs(i) += eps;
297
298 Eigen::VectorXd f2 = expansion2->Evaluate(evalPt, newCoeffs).at(0);
299
300 for(int j=0; j<f2.size(); ++j)
301 EXPECT_NEAR((f2(j)-f1(j))/eps, jac(j,i), 1e-3);
302 }
303 }
304
TEST_F(Approximation_MonotoneExpansion2d,Determinant)305 TEST_F(Approximation_MonotoneExpansion2d, Determinant){
306
307 Eigen::VectorXd evalPt(2);
308 evalPt << 0.1, 0.2;
309
310 Eigen::VectorXd coeffs = expansion1->GetCoeffs();
311
312 double logDet1 = expansion1->LogDeterminant(evalPt);
313
314 Eigen::VectorXd grad = expansion1->GradLogDeterminant(evalPt);
315
316 const double eps = 1e-7;
317
318 for(int i=0; i<coeffs.size(); ++i){
319 Eigen::VectorXd newCoeffs = coeffs;
320 newCoeffs(i) += eps;
321 double logDet2 = expansion1->LogDeterminant(evalPt,newCoeffs);
322
323 EXPECT_NEAR((logDet2-logDet1)/eps, grad(i), 1e-3);
324 }
325
326 }
327