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