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