1 #include "LegendreBasisForMinimumFrobeniusNormModel.hpp"
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 
6 
7 //--------------------------------------------------------------------------------
LegendreBasisForMinimumFrobeniusNormModel(int dim_input)8 LegendreBasisForMinimumFrobeniusNormModel::LegendreBasisForMinimumFrobeniusNormModel
9   ( int dim_input ) :
10   BasisForSurrogateModelBaseClass ( dim_input ),
11   QuadraticLegendre ( dim_input )
12 {
13   nb_basis_functions = ( dim_input * ( dim_input + 3 ) + 2 ) / 2;
14   nb_nodes = 0;
15 }
16 //--------------------------------------------------------------------------------
17 
18 //--------------------------------------------------------------------------------
set_nb_nodes(int nb_nodes_input)19 void LegendreBasisForMinimumFrobeniusNormModel::set_nb_nodes( int nb_nodes_input ) {
20   if ( nb_nodes == nb_nodes_input ) return;
21   nb_nodes = nb_nodes_input;
22 
23   basis_values.clear();
24   basis_values.resize ( nb_nodes );
25   A_sysmat = Eigen::MatrixXd::Zero(nb_basis_functions + nb_nodes,
26                                    nb_basis_functions + nb_nodes);
27 
28   for ( int i = 0; i < nb_basis_functions; ++i)
29     if (i > BasisForSurrogateModelBaseClass::dim) A_sysmat(i, i) = 1e0;
30   F_rhsmat = Eigen::MatrixXd::Zero(nb_basis_functions + nb_nodes, nb_nodes);
31 
32   for ( int i = 0; i < nb_nodes; ++i )
33     F_rhsmat(nb_basis_functions + i, i) = 1e0;
34   basis_constants.clear();
35   basis_constants.resize( nb_nodes );
36   basis_gradients.clear();
37   basis_gradients.resize( nb_nodes );
38   basis_Hessians.clear();
39   basis_Hessians.resize( nb_nodes );
40   basis_coefficients.clear();
41   basis_coefficients.resize( nb_nodes );
42 
43 
44   for ( int i = 0; i < nb_nodes; ++i ) {
45     basis_coefficients[i] = Eigen::VectorXd::Zero( nb_basis_functions ) ;
46     basis_gradients[i].resize( BasisForSurrogateModelBaseClass::dim );
47     basis_Hessians[i].resize( BasisForSurrogateModelBaseClass::dim );
48     for ( int j = 0; j < BasisForSurrogateModelBaseClass::dim; ++j)
49       basis_Hessians[i][j].resize( BasisForSurrogateModelBaseClass::dim );
50   }
51 
52   return;
53 }
54 //--------------------------------------------------------------------------------
55 
56 
57 //--------------------------------------------------------------------------------
compute_basis_coefficients(std::vector<std::vector<double>> const & nodes)58 void LegendreBasisForMinimumFrobeniusNormModel::compute_basis_coefficients (
59   std::vector< std::vector<double> > const &nodes )
60 {
61 
62   set_nb_nodes ( nodes.size( ) );
63 
64 
65   // system matrix for computing coeffs of Lagrange interpolation models
66   for ( int i = 0; i < nb_basis_functions; ++i ) {
67     for ( int j = 0; j < nb_nodes; ++j ) {
68       A_sysmat(i, j+nb_basis_functions) = evaluate_basis( i, nodes[j] );
69       A_sysmat(j+nb_basis_functions, i) = A_sysmat(i, j+nb_basis_functions);
70     }
71   }
72 
73 
74   S_coeffsolve = A_sysmat.colPivHouseholderQr().solve(F_rhsmat);
75 
76   if ( (A_sysmat * S_coeffsolve - F_rhsmat).norm() > 1e-5) {
77     for ( int i = 0; i < nb_nodes; ++i ) {
78       for ( int j = 0; j < BasisForSurrogateModelBaseClass::dim; ++j )
79         S_coeffsolve(BasisForSurrogateModelBaseClass::dim +j+1, i) = 0e0;
80       if ( norm( nodes[i]) < 1e-16 )
81         S_coeffsolve(0, i) = 1e0;
82       else
83         S_coeffsolve(0, i) = 0e0;
84     }
85   }
86 
87   for ( int i = 0; i < nb_nodes; ++i ) {
88     basis_coefficients[i] = S_coeffsolve.block(0,0,nb_basis_functions, nb_nodes).col(i);
89     compute_mat_vec_representation ( i );
90   }
91 
92 
93 
94 
95 
96   return;
97 }
98 //--------------------------------------------------------------------------------
99 
100 
101 //--------------------------------------------------------------------------------
value(int basis_number)102 double &LegendreBasisForMinimumFrobeniusNormModel::value ( int basis_number )
103 {
104   return basis_constants[ basis_number ];
105 }
106 //--------------------------------------------------------------------------------
107 
108 //--------------------------------------------------------------------------------
gradient(int basis_number)109 std::vector<double> &LegendreBasisForMinimumFrobeniusNormModel::gradient (int basis_number)
110 {
111   return basis_gradients.at( basis_number );
112 }
113 //--------------------------------------------------------------------------------
114 
115 //--------------------------------------------------------------------------------
hessian(int basis_number)116 std::vector< std::vector<double> > &LegendreBasisForMinimumFrobeniusNormModel::hessian (
117   int basis_number )
118 {
119   return basis_Hessians[ basis_number ];
120 }
121 //--------------------------------------------------------------------------------
122 
123 
124 //--------------------------------------------------------------------------------
compute_mat_vec_representation(int basis_number)125 void LegendreBasisForMinimumFrobeniusNormModel::compute_mat_vec_representation ( int basis_number )
126 // computes the representation m(x) = c + g.dot(x) + 0.5*x.dot(H*x) in scaled form
127 // x_best = 0 and delta = 1;
128 {
129   counter = 0;
130   basis_constants[basis_number] = basis_coefficients[basis_number](0);
131   for (int j = 0; j < BasisForSurrogateModelBaseClass::dim; ++j) {
132     basis_constants[basis_number] -= 0.5*basis_coefficients[basis_number](
133                                          j+1+BasisForSurrogateModelBaseClass::dim);
134     basis_gradients[basis_number].at(j) = basis_coefficients[ basis_number ]( j + 1 );
135     basis_Hessians[basis_number].at(j).at(j) = 3e0*basis_coefficients[ basis_number ]
136                                                (j+1+BasisForSurrogateModelBaseClass::dim);
137     for (int k = j+1; k < BasisForSurrogateModelBaseClass::dim; ++k) {
138       basis_Hessians[basis_number].at(j).at(k) = basis_coefficients[ basis_number ]
139                                                  (2*BasisForSurrogateModelBaseClass::dim+1+counter);
140       basis_Hessians[basis_number].at(k).at(j) = basis_coefficients[ basis_number ]
141                                                  (2*BasisForSurrogateModelBaseClass::dim+1+counter);
142       counter++;
143     }
144   }
145 
146   return;
147 }
148 //--------------------------------------------------------------------------------
149 
150 
151 //--------------------------------------------------------------------------------
evaluate(std::vector<double> const & x)152 std::vector<double> &LegendreBasisForMinimumFrobeniusNormModel::evaluate (
153   std::vector<double> const &x )
154 {
155 
156   //assert( nb_nodes == basis_coefficients.size());
157 
158   for ( int i = 0; i < nb_nodes; ++i ) {
159     basis_values.at( i ) = 0e0;
160     for ( int j = 0; j < nb_basis_functions; ++j ) {
161       basis_values.at( i ) += basis_coefficients[ i ]( j ) * evaluate_basis ( j, x );
162     }
163   }
164 
165   return basis_values;
166 }
167 //--------------------------------------------------------------------------------
168 
169 
170 //--------------------------------------------------------------------------------
evaluate(std::vector<double> const & x,int basis_number)171 double LegendreBasisForMinimumFrobeniusNormModel::evaluate (
172   std::vector<double> const &x, int basis_number)
173 {
174 
175   double basis_value = 0e0;
176   for ( int j = 0; j < nb_basis_functions; ++j ) {
177     basis_value += basis_coefficients[ basis_number ]( j ) *
178                    evaluate_basis ( j, x );
179   }
180   return basis_value;
181 }
182 //--------------------------------------------------------------------------------
183 
184 
185