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