1 #ifndef HPostProcessModels 2 #define HPostProcessModels 3 4 #include "BlackBoxData.hpp" 5 #include "MinimumFrobeniusNormModel.hpp" 6 #include "LegendreBasisForMinimumFrobeniusNormModel.hpp" 7 #include "VectorOperations.hpp" 8 9 template<class TSurrogateModel = MinimumFrobeniusNormModel, 10 class TBasisForSurrogateModel = LegendreBasisForMinimumFrobeniusNormModel> 11 class PostProcessModels : VectorOperations { 12 private: 13 int nb_models; 14 int dim; 15 double delta; 16 std::unique_ptr<TSurrogateModel> surrogate_model_prototype; 17 TBasisForSurrogateModel surrogate_basis; 18 std::vector<TSurrogateModel> surrogate_models; 19 std::shared_ptr<BlackBoxData> evaluations; 20 std::vector<double> vec_output; 21 std::vector<double> c_vec; 22 std::vector< std::vector<double> > g_vec; 23 std::vector< std::vector< std::vector<double> > > H_vec; 24 public: PostProcessModels(BlackBoxData & evaluations_input)25 PostProcessModels ( BlackBoxData &evaluations_input ) : 26 surrogate_basis ( evaluations_input.nodes[0].size() ) 27 { 28 surrogate_model_prototype.reset( new TSurrogateModel(surrogate_basis) ); 29 evaluations = std::make_shared<BlackBoxData>(evaluations_input); 30 nb_models = evaluations->values.size(); 31 dim = evaluations->nodes[0].size(); 32 delta = evaluations->get_scaling(); 33 vec_output.resize(dim); 34 c_vec.resize( nb_models ); 35 g_vec.resize( nb_models ); 36 H_vec.resize( nb_models ); 37 38 surrogate_basis.compute_basis_coefficients ( evaluations->get_scaled_active_nodes( delta ) ) ; 39 40 for (int i = 0; i < nb_models; ++i ) { 41 surrogate_models.push_back ( *surrogate_model_prototype ); 42 surrogate_models[ i ].set_function_values ( evaluations->get_active_values(i) ); 43 } 44 45 for ( int i = 0; i < nb_models; ++i ) { 46 g_vec[i].resize( dim ); 47 H_vec[i].resize( dim ); 48 for ( int j = 0; j < dim; ++j ) { 49 H_vec[i][j].resize( dim ); 50 } 51 get_c_g_H( i, c_vec[i], g_vec[i], H_vec[i]); 52 } 53 54 return; 55 } 56 57 //-------------------------------------------------------------------------------- get_c_g_H(int model_number,double & c,std::vector<double> & g,std::vector<std::vector<double>> & H)58 void get_c_g_H ( int model_number, double &c, std::vector<double> &g, 59 std::vector< std::vector<double> > &H ) 60 { 61 62 double evaluations_scaling = delta; 63 std::vector<double> evaluations_center_node = evaluations->nodes[ evaluations->best_index ]; 64 for (int i = 0; i < dim; ++i) g[i] = 0.0; 65 c = surrogate_models[model_number].evaluate( g ); 66 c -= VectorOperations::dot_product( surrogate_models[model_number].gradient(), 67 evaluations_center_node) / evaluations_scaling; 68 this->mat_vec_product( surrogate_models[model_number].hessian(), 69 evaluations_center_node, g ); 70 c += 0.5*VectorOperations::dot_product( evaluations_center_node, g ) / 71 ( evaluations_scaling * evaluations_scaling ); 72 73 this->scale ( -1.0/(evaluations_scaling * evaluations_scaling), g, g); 74 this->add( 1.0/evaluations_scaling, surrogate_models[model_number].gradient(), g); 75 76 H = surrogate_models[model_number].hessian(); 77 for ( int i = 0; i < dim; ++i ) { 78 for ( int j = 0; j < dim; ++j ) { 79 H[i][j] /= (evaluations_scaling * evaluations_scaling); 80 } 81 } 82 83 return; 84 } 85 //-------------------------------------------------------------------------------- 86 //-------------------------------------------------------------------------------- get_trustregion()87 double get_trustregion ( ) 88 { 89 return evaluations->get_scaling(); 90 } 91 //-------------------------------------------------------------------------------- 92 //-------------------------------------------------------------------------------- get_c(int model_number,std::vector<double> & x)93 double get_c ( int model_number, std::vector<double> &x ) 94 { 95 assert( x.size() == (unsigned) dim ); 96 double c = c_vec[ model_number ]; 97 c += VectorOperations::dot_product( g_vec[ model_number], x ); 98 this->mat_vec_product( H_vec[model_number], x, vec_output ); 99 c += 0.5*VectorOperations::dot_product( vec_output, x ); 100 return c; 101 } 102 //-------------------------------------------------------------------------------- 103 //-------------------------------------------------------------------------------- get_g(int model_number,std::vector<double> & x)104 std::vector<double> get_g ( int model_number, std::vector<double> &x ) 105 { 106 this->mat_vec_product( H_vec[model_number], x, vec_output); 107 this->add( g_vec[model_number], vec_output ); 108 return vec_output; 109 } 110 //-------------------------------------------------------------------------------- 111 //-------------------------------------------------------------------------------- get_H(int model_number)112 std::vector< std::vector<double> > get_H ( int model_number ) 113 { 114 return H_vec[model_number]; 115 } 116 //-------------------------------------------------------------------------------- 117 118 }; 119 120 #endif 121