1 /*-------------------------------------------------------------------------------------*/
2 /*  NOMAD - Nonlinear Optimization by Mesh Adaptive Direct search - version 3.7.2      */
3 /*                                                                                     */
4 /*  Copyright (C) 2001-2015  Mark Abramson        - the Boeing Company, Seattle        */
5 /*                           Charles Audet        - Ecole Polytechnique, Montreal      */
6 /*                           Gilles Couture       - Ecole Polytechnique, Montreal      */
7 /*                           John Dennis          - Rice University, Houston           */
8 /*                           Sebastien Le Digabel - Ecole Polytechnique, Montreal      */
9 /*                           Christophe Tribes    - Ecole Polytechnique, Montreal      */
10 /*                                                                                     */
11 /*  funded in part by AFOSR and Exxon Mobil                                            */
12 /*                                                                                     */
13 /*  Author: Sebastien Le Digabel                                                       */
14 /*                                                                                     */
15 /*  Contact information:                                                               */
16 /*    Ecole Polytechnique de Montreal - GERAD                                          */
17 /*    C.P. 6079, Succ. Centre-ville, Montreal (Quebec) H3C 3A7 Canada                  */
18 /*    e-mail: nomad@gerad.ca                                                           */
19 /*    phone : 1-514-340-6053 #6928                                                     */
20 /*    fax   : 1-514-340-5665                                                           */
21 /*                                                                                     */
22 /*  This program is free software: you can redistribute it and/or modify it under the  */
23 /*  terms of the GNU Lesser General Public License as published by the Free Software   */
24 /*  Foundation, either version 3 of the License, or (at your option) any later         */
25 /*  version.                                                                           */
26 /*                                                                                     */
27 /*  This program is distributed in the hope that it will be useful, but WITHOUT ANY    */
28 /*  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A    */
29 /*  PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more details.   */
30 /*                                                                                     */
31 /*  You should have received a copy of the GNU Lesser General Public License along     */
32 /*  with this program. If not, see <http://www.gnu.org/licenses/>.                     */
33 /*                                                                                     */
34 /*  You can find information on the NOMAD software at www.gerad.ca/nomad               */
35 /*-------------------------------------------------------------------------------------*/
36 /**
37   \file   Quad_Model.hpp
38   \brief  Quadratic regression or MFN interpolation model (headers)
39   \author Sebastien Le Digabel
40   \date   2010-08-31
41   \see    Quad_Model.cpp
42 */
43 #ifndef __QUAD_MODEL__
44 #define __QUAD_MODEL__
45 
46 #include "Cache.hpp"
47 #include "Model_Sorted_Point.hpp"
48 #include "Evaluator.hpp"
49 
50 namespace NOMAD {
51 
52   /// Class for quadratic regression or MFN interpolation model.
53   class Quad_Model : private NOMAD::Uncopyable {
54 
55     /*-------------------------------------------------------------------------*/
56   private:
57 
58     const NOMAD::Display                     & _out;  ///< Display.
59 
60     std::vector<NOMAD::Eval_Point *>           _Y;    ///< Interpolation points.
61 
62     const std::vector<NOMAD::bb_output_type> & _bbot; ///< Blackbox output types.
63 
64     NOMAD::interpolation_type _interpolation_type;    ///< Interpolation type.
65 
66     int                       _n;                     ///< Dimension.
67     int                       _nfree;                 ///< Number of free variables.
68     int                       _n_alpha;               ///< Number of model coefficients.
69     bool                    * _fixed_vars;            ///< Fixed variables.
70     int                     * _index;                 ///< Var. indexes with fixed var.
71     NOMAD::Point           ** _alpha;                 ///< Model coefficients.
72     NOMAD::Point              _center;                ///< Model center.
73     NOMAD::Point              _ref;                   ///< Reference for scaling.
74     NOMAD::Point              _scaling;               ///< Scaling.
75     const NOMAD::Cache      & _cache;                 ///< Cache.
76     const NOMAD::Signature  & _signature;             ///< Signature.
77     bool                      _error_flag;            ///< Error flag.
78 	std::list<NOMAD::Direction>  _dirP;               ///< Directions used for scaling (may be empty)
79 	NOMAD::Point _delta_m;                           ///< Mesh size used for scaling
80 	NOMAD::Double _epsilon;                           ///< Offset for direction scaling
81 
82 
83     NOMAD::Double             _cond;                  ///< Condition number.
84 
85     /// Initialize alpha (model parameters).
86     void init_alpha ( void );
87 
88     /// Check if an unscaled point is in \c B(center,radius) for a given radius.
89     /**
90        \param x      The unscaled point -- \b IN.
91        \param radius The radius         -- \b IN.
92        \return \c true is \c x is in \c B(center,radius).
93     */
94     bool is_within_radius ( const NOMAD::Point & x      ,
95 			    const NOMAD::Point & radius   ) const;
96 
97     /// Check the interpolation set \c Y.
98     /**
99        \return \c true if the interpolation set is valid.
100     */
101     bool check_Y ( void ) const;
102 
103     /// Check outputs before the integration into \c Y.
104     /**
105        \param bbo The outputs       -- \b IN.
106        \param m   Number of outputs -- \b IN.
107        return \c true if the \c m outputs are valid.
108     */
109     bool check_outputs ( const NOMAD::Point & bbo , int m ) const;
110 
111     /// Reduce the number of interpolation points.
112     /**
113        The points are sorted accorded to their distance to the model center.
114        \param center      Center of the model                -- \b IN.
115        \param max_Y_size  Max number of interpolation points -- \b IN.
116     */
117     void reduce_Y  ( const NOMAD::Point & center , int max_Y_size );
118 
119     /// Compute condition number.
120     /**
121        \param W   Matrix W given as a vector -- \b IN.
122        \param n   Size of \c W               -- \b IN
123        \param eps Epsilon                    -- \b IN.
124     */
125     void compute_cond ( const double * W , int n , double eps );
126 
127     /// Compute the cumulated error of a model for one output.
128     /**
129        The errors are computed on the interpolation set \c Y.
130        \param bbo_index   Blackbox output index  -- \b IN.
131        \param error       Cumulated error        -- \b OUT.
132        \param min_rel_err Min relative error     -- \b OUT.
133        \param max_rel_err Max relative error     -- \b OUT.
134        \param avg_rel_err Average relative error -- \b OUT.
135     */
136     void compute_model_error ( int             bbo_index   ,
137 			       NOMAD::Double & error       ,
138 			       NOMAD::Double & min_rel_err ,
139 			       NOMAD::Double & max_rel_err ,
140 			       NOMAD::Double & avg_rel_err   ) const;
141 
142     /// Compute the maximal relative error of a model.
143     /**
144        The error is computed for the interpolation set \c Y.
145        \return The maximal relative error.
146     */
147     NOMAD::Double compute_max_rel_err ( void ) const;
148 
149     /// Compute the element \c (i,j) of the interpolation matrix \c M(phi,Y).
150     /**
151        \param i Row index    -- \b IN.
152        \param j Column index -- \b IN.
153     */
154     double compute_M ( int i , int j ) const;
155 
156     /// Construct Minimum Frobenius Norm (MFN) model.
157     /**
158        - This occurs when \c p+1 \c < \c (n+1)(n+2)/2.
159        \param eps        Epsilon                               -- \b IN.
160        \param max_mpn    Maximum \c m+n value for SVD matrices -- \b IN.
161        \param max_Y_size Maximum number of elements in \c Y    -- \b IN.
162        \return \c true if the construction succeeded
163     */
164     bool construct_MFN_model ( double eps , int max_mpn , int max_Y_size );
165 
166     /// Construct regression model.
167     /**
168        - This occurs when \c p+1 \c >= \c (n+1)(n+2)/2.
169        \param eps        Epsilon                               -- \b IN.
170        \param max_mpn    Maximum \c m+n value for SVD matrices -- \b IN.
171        \param max_Y_size Maximum number of elements in \c Y    -- \b IN.
172        \return \c true if the construction succeeded
173     */
174     bool construct_regression_model ( double eps        ,
175 				      int    max_mpn    ,
176 				      int    max_Y_size   );
177 
178     /// Construct well-poised (WP) model.
179     /**
180        \param max_Y_size Maximum number of elements in \c Y -- \b IN.
181        \return \c true if the construction succeeded
182     */
183     bool construct_WP_model ( int max_Y_size );
184 
185     /// Find interpolation point with max Lagrange polynomial value.
186     /**
187        \param  li      Lagrange polynomial             -- \b IN.
188        \param  Y       Interpolation points candidates -- \b IN.
189        \param  i1      Initial index in \c Y           -- \b IN.
190        \param  i2      Final index in \c Y             -- \b IN.
191        \param  max_lix Absolute value of the max value -- \b OUT.
192        \return Index of interpolation point.
193     */
194     int find_max_lix ( const NOMAD::Point                     & li      ,
195 		       const std::vector<NOMAD::Eval_Point *> & Y       ,
196 		       int                                      i1      ,
197 		       int                                      i2      ,
198 		       NOMAD::Double                          & max_lix   ) const;
199 
200     /// Resolution of system \c F.[mu alpha_L]'=[f(Y) 0]' for MFN interpolation.
201     /**
202        \param U         Matrix \c F=U from the SVD decomposition \c U.W.V' -- \b IN.
203        \param W         Matrix \c W from the SVD decomposition \c U.W.V'   -- \b IN.
204        \param V         Matrix \c V from the SVD decomposition \c U.W.V'   -- \b IN.
205        \param bbo_index Blackbox output index                              -- \b IN.
206        \param alpha     Model parameters                                   -- \b IN.
207        \param eps       Epsilon                                            -- \b IN.
208     */
209     void solve_MFN_system ( double      ** U         ,
210 			    double       * W         ,
211 			    double      ** V         ,
212 			    int            bbo_index ,
213 			    NOMAD::Point & alpha     ,
214 			    double         eps	      ) const;
215 
216     /// Resolution of system \c F.alpha=M'.f(Y) for the regression.
217     /**
218        \param M         Matrix \c M                                        -- \b IN.
219        \param U         Matrix \c F=U from the SVD decomposition \c U.W.V' -- \b IN.
220        \param W         Matrix \c W from the SVD decomposition \c U.W.V'   -- \b IN.
221        \param V         Matrix \c V from the SVD decomposition \c U.W.V'   -- \b IN.
222        \param bbo_index Blackbox output index                              -- \b IN.
223        \param alpha     Model parameters                                   -- \b IN.
224        \param eps       Epsilon                                            -- \b IN.
225     */
226     void solve_regression_system ( double      ** M         ,
227 				   double      ** U         ,
228 				   double       * W         ,
229 				   double      ** V         ,
230 				   int            bbo_index ,
231 				   NOMAD::Point & alpha     ,
232 				   double         eps	      ) const;
233 
234     /// Display Lagrange polynomials.
235     /**
236        \param l Lagrange polynomials -- \b IN.
237        \param Y Interpolation set    -- \b IN.
238     */
239     void display_lagrange_polynomials
240     ( const std::vector<NOMAD::Point      *> & l ,
241       const std::vector<NOMAD::Eval_Point *> & Y   ) const;
242 
243 #ifdef MODEL_STATS
244     mutable NOMAD::Double _Yw; ///< Width of the interpolation set \c Y.
245 
246   public:
247 
248     /// Access to the width of the interpolation set \c X (or \c Y).
249     /**
250        \return The width of the interpolation set.
251     */
get_Yw(void) const252     const NOMAD::Double & get_Yw ( void ) const { return _Yw; }
253 #endif
254 
255     /*-------------------------------------------------------------------------*/
256   public:
257 
258     /// Constructor.
259     /**
260        \param out           The NOMAD::Display object   -- \b IN.
261        \param bbot          Output types                -- \b IN.
262        \param cache         Cache                       -- \b IN.
263        \param signature     Signature                   -- \b IN.
264     */
265     Quad_Model ( const NOMAD::Display                     & out       ,
266 		 const std::vector<NOMAD::bb_output_type> & bbot      ,
267 		 const NOMAD::Cache                       & cache     ,
268 		 const NOMAD::Signature                   & signature   );
269 
270     /// Destructor.
271     virtual ~Quad_Model ( void );
272 
273     /// Evaluate a model at a given point.
274     /**
275        \param x     The point        -- \b IN.
276        \param alpha Model parameters -- \b IN.
277        \return Model value.
278     */
279     NOMAD::Double eval ( const NOMAD::Point & x     ,
280 			 const NOMAD::Point & alpha   ) const;
281 
282 	  /// Compute model \c h and \c f values at a point.
283     /**
284        \param x      The point                 -- \b IN.
285        \param h_min  Value of \c h_min         -- \b IN..
286        \param h_norm Norm used to compute \c h -- \b IN..
287        \param h      Value of \c h             -- \b OUT.
288        \param f      Value of \c f             -- \b OUT.
289     */
290     void eval_hf ( const NOMAD::Point  & x      ,
291 		   const NOMAD::Double & h_min  ,
292 		   NOMAD::hnorm_type     h_norm ,
293 		   NOMAD::Double       & h      ,
294 		   NOMAD::Double       & f        ) const;
295 
296 
297     /// Access to the interpolation type.
298     /**
299        \return The interpolation type.
300     */
get_interpolation_type(void) const301     const NOMAD::interpolation_type & get_interpolation_type ( void ) const
302     {
303       return _interpolation_type;
304     }
305 
306     /// Access to the center of the model.
307     /**
308        \return The center.
309     */
get_center(void) const310     const NOMAD::Point & get_center ( void ) const { return _center; }
311 
312     /// Access to the dimension.
313     /**
314        \return The dimension \c n.
315     */
get_n(void) const316     int get_n ( void ) const { return _n; }
317 
318     /// Access to the number of free variables.
319     /**
320        \return The number of free variables \c n.
321     */
get_nfree(void) const322     int get_nfree ( void ) const { return _nfree; }
323 
324     /// Access to the model parameters.
325     /**
326        \return The model parameters \c alpha.
327     */
get_alpha(void) const328     NOMAD::Point ** get_alpha ( void ) const { return _alpha; }
329 
330     /// Check if the model is ready for evaluations.
331     /**
332        \return A boolean equal to \c true if the model is ready.
333     */
334     bool check ( void ) const;
335 
336     /// Access to the fixed variables.
337     /**
338        \param i Variable index -- \b IN.
339        \return \c true if variable \c i is fixed.
340     */
variable_is_fixed(int i) const341     bool variable_is_fixed ( int i ) const { return _fixed_vars[i]; }
342 
343     /// Access to the number of interpolation points.
344     /**
345        \return The number of interpolation points \c nY=p+1.
346     */
get_nY(void) const347     int get_nY ( void ) const { return static_cast<int> ( _Y.size() ); }
348 
349     /// Access to the condition number.
350     /**
351        \return The condition number.
352     */
get_cond(void) const353     const NOMAD::Double & get_cond ( void ) const { return _cond; }
354 
355     /// Access to the error flag.
356     /**
357        \return The error flag.
358     */
get_error_flag(void) const359     bool get_error_flag ( void ) const { return _error_flag; }
360 
361     /// Construct the interpolation set \c Y.
362     /**
363        \param center               Model center                       -- \b IN.
364        \param interpolation_radius Interpolation radius               -- \b IN.
365        \param max_Y_size           Maximum number of elements in \c Y -- \b IN.
366     */
367     void construct_Y ( const NOMAD::Point & center               ,
368 		       const NOMAD::Point & interpolation_radius ,
369 		       int                  max_Y_size             );
370 
371     /// Construct \c m models (one by output).
372     /**
373        \param use_WP     Use or not well-poisedness            -- \b IN.
374        \param eps        Epsilon                               -- \b IN.
375        \param max_mpn    Maximum \c m+n value for SVD matrices -- \b IN.
376        \param max_Y_size Maximum number of elements in \c Y    -- \b IN.
377     */
378     void construct ( bool   use_WP     ,
379 		     double eps        ,
380 		     int    max_mpn    ,
381 		     int    max_Y_size   );
382 
383     /// Define scaling to put all coordinates centered in \c [-r;r].
384     /**
385        - Looks also for fixed variables.
386        \param r The \c r parameter corresponds to \c MODEL_RADIUS_FACTOR -- \b IN.
387     */
388     void define_scaling ( const NOMAD::Double & r );
389 
390 	  /// Define scaling based on directions. See paper: Reducing the number of function evaluations in Mesh Adaptive Direct Search algorithms, Audet, Ianni, LeDigabel, Tribes, 2014
391 	  /**
392        - Looks also for fixed variables.
393        \param dirP    The \c dirP parameter corresponds to set of directions formin a hyper-cube centered on poll center -- \b IN.
394 	   \param delta_m The \c delta_m parameter is the dimension of the mesh -- \b IN.
395 	   \param epsilon The \c epsilon parameter is the hyper-cube offset from the poll center -- \b IN.
396 	   */
397 	void define_scaling_by_directions ( const std::list<NOMAD::Direction> & dirP, const NOMAD::Point & delta_m, const NOMAD::Double &epsilon  );
398 
399 
400     /// Scale a point.
401     /**
402        \param x The point to scale -- \b IN/OUT.
403        \return \c true if the scaling worked.
404     */
405     bool scale ( NOMAD::Point & x ) const;
406 
407     /// Unscale a point.
408     /**
409        \param x The point to unscale -- \b IN/OUT.
410        \return \c true if the unscaling worked.
411     */
412     bool unscale ( NOMAD::Point & x ) const;
413 
414 	  /// Unscale the gradient at a point.
415 	  /**
416        \param x The grad to unscale -- \b IN/OUT.
417        \return \c true if the unscaling worked.
418 	   */
419 	  bool unscale_grad ( NOMAD::Point & x ) const;
420 
421 
422     /// Check if a caled point is inside the trust radius.
423     /**
424        \param x The scaled point -- \b IN.
425        \return  \c true is \c x is in the trust radius.
426     */
427     bool is_within_trust_radius ( const NOMAD::Point & x ) const;
428 
429     /// Display the model coefficients.
430     /**
431        \param out   The NOMAD::Display object -- \b IN.
432     */
433     void display_model_coeffs ( const NOMAD::Display & out ) const;
434 
435     /// Display the interpolation set \c Y.
436     /**
437        \param out   The NOMAD::Display object  -- \b IN.
438        \param title Title of the display block -- \b IN
439                     --\b optional (default="interpolation set Y").
440     */
441     void display_Y ( const NOMAD::Display & out ,
442 		     const std::string    & title = "interpolation set Y" ) const;
443 
444     /// Display cumulated error on the interpolation points.
445     /**
446        \param out   The NOMAD::Display object -- \b IN.
447     */
448     void display_Y_error ( const NOMAD::Display & out ) const;
449   };
450 }
451 
452 #endif
453