/*-------------------------------------------------------------------------------------*/
/*  NOMAD - Nonsmooth Optimization by Mesh Adaptive Direct search - version 3.0.0      */
/*                                                                                     */
/*  Copyright (C) 2001-2008  Mark Abramson        - the Boeing Company, Seattle        */
/*                           Charles Audet        - Ecole Polytechnique, Montreal      */
/*                           Gilles Couture       - Ecole Polytechnique, Montreal      */
/*                           John Dennis          - Rice University, Houston           */
/*                           Sebastien Le Digabel - Ecole Polytechnique, Montreal      */
/*                                                                                     */
/*  funded in part by AFOSR and Exxon Mobil                                            */
/*                                                                                     */
/*  Author: Sebastien Le Digabel                                                       */
/*                                                                                     */
/*  Contact information:                                                               */
/*    Ecole Polytechnique de Montreal - GERAD                                          */
/*    C.P. 6079, Succ. Centre-ville, Montreal (Quebec) H3C 3A7 Canada                  */
/*    e-mail: Sebastien.Le.Digabel@gerad.ca                                            */
/*    phone : 1-514-340-6053 #6673                                                     */
/*    fax   : 1-514-340-5665                                                           */
/*                                                                                     */
/*  This program is free software: you can redistribute it and/or modify it under the  */
/*  terms of the GNU General Public License as published by the Free Software          */
/*  Foundation, either version 3 of the License, or (at your option) any later         */
/*  version.                                                                           */
/*                                                                                     */
/*  This program is distributed in the hope that it will be useful, but WITHOUT ANY    */
/*  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A    */
/*  PARTICULAR PURPOSE.  See the GNU General Public License for more details.          */
/*                                                                                     */
/*  You should have received a copy of the GNU General Public License along with this  */
/*  program. If not, see <http://www.gnu.org/licenses/>.                               */
/*                                                                                     */
/*  You can find information on the NOMAD software at www.gerad.ca/nomad               */
/*-------------------------------------------------------------------------------------*/
#include "Mads.hpp"
using namespace std;

/*-----------------------------------*/
/*   static members initialization   */
/*-----------------------------------*/
bool Mads::_force_quit = false;

/*---------------------------------------------------------*/
/*                         constructor                     */
/*---------------------------------------------------------*/
Mads::Mads ( Parameters & p , Evaluator * ev , ostream & out ) :
  _p             ( p                     ) ,
  _stats         ( p.get_sgte_cost()     ) ,
  _std_signature ( p , out               ) ,
  _ev_control    ( p , _stats , ev , out ) ,
  _barrier       ( p , out               ) ,
  _out           ( out                   ) {

  // random seed initialization:
  srand ( p.get_seed() );
}

/*---------------------------------------------------------*/
/*           force quit (called by pressing ctrl-c)        */
/*---------------------------------------------------------*/
void Mads::force_quit ( int signalValue ) {
  if ( signalValue == SIGINT ) {
    Mads::_force_quit = true;
    Evaluator_Control::force_quit();
  }
}

/*---------------------------------------------------------*/
/*                         reset                           */
/*---------------------------------------------------------*/
void Mads::reset ( bool keep_barrier , bool keep_stats ) {

  // barrier:
  if ( !keep_barrier )
    _barrier.reset();

  // standard signature:
  _std_signature.reset(_p);

  // stats:
  if ( !keep_stats )
    _stats.reset();

  // mesh:
  Mesh::init ( _p.get_mesh_update_basis().value() ,
	       _p.get_mesh_coarsening_exponent () ,
	       _p.get_mesh_refining_exponent   () ,
	       _p.get_initial_mesh_index       ()   );
}

/*---------------------------------------------------------*/
/*                    algorithm execution                  */
/*---------------------------------------------------------*/
stop_type Mads::run ( void ) {

  // displays:
  int display_degree = _p.get_display_degree();
  if ( display_degree > 2 )
    _out << "\n\nbegin of Mads::run()\n--------------------\n\n";
  else if ( display_degree == 2 ) {
    _out << "\nbegin of Mads::run()\n\n";
    _ev_control.display_stats ( true                   ,
				_out                   ,
				_p.get_display_stats() ,
				NOMAD::Double()        ,
				NOMAD::Double()        ,
				NULL                   ,
				NULL                     );
    _out << endl << endl;
  }

  // barrier init:
  _barrier.reset();

  // mesh init:
  Mesh::init ( _p.get_mesh_update_basis().value() ,
	       _p.get_mesh_coarsening_exponent () ,
	       _p.get_mesh_refining_exponent   () ,
	       _p.get_initial_mesh_index       ()   );

  bool               stop               = false;
  stop_type          stop_reason        = _NO_STOP_;
  const Eval_Point * new_feas_inc       = NULL;
  const Eval_Point * new_infeas_inc     = NULL;

  Direction  feas_success_dir;
  Direction  infeas_success_dir;

  // x0 eval:
  eval_x0 ( stop , stop_reason );

  // phase one: if no feasible starting point:
  if ( stop                     &&
       stop_reason == _X0_EVAL_ &&
       _stats.get_eval() > 0    &&
       !phase_one ( stop , stop_reason ) ) {
    stop        = true;
    stop_reason = _PHASE_1_FAIL_;
  }

  // initial Latin-Hypercube (LH) search:
  if ( _p.get_LH_search_p0() > 0 ) {
    success_type success;
    LH_search ( true           ,
		_std_signature ,
		stop           ,
		stop_reason    ,
		success        ,
		new_feas_inc   ,
		new_infeas_inc   );
  }

  // no iterations allowed:
  if ( !stop && _p.get_max_iterations() == 0 ) {
    stop        = true;
    stop_reason = _MAX_ITER_REACHED_;
  }

  // MADS iterations:
  while ( !stop )
    iteration ( stop               ,
		stop_reason        ,
		feas_success_dir   ,
		infeas_success_dir ,
		new_feas_inc       ,
		new_infeas_inc       );

  // final cache save (overwrite=true):
  _ev_control.save_caches(true);

  // final displays:
  const Eval_Point * bf = _barrier.get_best_feasible();

  const string stats_file_name = _p.get_stats_file_name();
  if ( bf && !stats_file_name.empty() )
    _ev_control.stats_file ( stats_file_name       ,
			     bf->get_f()           ,
			     bf->get_h()           ,
			     bf                    ,
			     &bf->get_bb_outputs()   );

  if ( display_degree > 0 ) {

    // final stats:
    if ( display_degree < 3 ) {

      if ( bf )
	_ev_control.display_stats ( false                  ,
				    _out                   ,
				    _p.get_display_stats() ,
				    bf->get_f()            ,
				    bf->get_h()            ,
				    bf                     ,
				    &bf->get_bb_outputs()    );
      else if ( display_degree > 1 ) {
	_out << "no feasible solution has been found" << endl;
	bf = _barrier.get_best_infeasible();
	if ( bf ) {
	  _out << "best infeasible solution (h=" << bf->get_h() << "):\n\t";
	  _ev_control.display_stats ( false                  ,
				      _out                   ,
				      _p.get_display_stats() ,
				      bf->get_f()            ,
				      bf->get_h()            ,
				      bf                     ,
				      &bf->get_bb_outputs()    );
	}
      }
    }

    if ( display_degree > 1 )
      _out << endl << "end of Mads::run(), " << stop_reason << endl;
  }

  return stop_reason;
}

/*---------------------------------------------------------*/
/*        MADS phase 1 (try to satisfy EB constraints)     */
/*---------------------------------------------------------*/
bool Mads::phase_one ( bool & stop , stop_type & stop_reason ) {

  int display_degree = _p.get_display_degree();
  if ( display_degree > 2 )
    _out << endl << "begin of Mads phase one" << endl;

  // counters:
  int old_bbe = _stats.get_bb_eval();
  int old_it  = _stats.get_iterations();

  // modify parameters:
  const string      old_stats_file_name = _p.get_stats_file_name();
  const list<string>     old_stats_file = _p.get_stats_file();
  const NOMAD::Double      old_f_target = _p.get_f_target();
  const vector<bb_output_type> old_bbot = _p.get_bb_output_types();
  vector<bb_output_type>        p1_bbot = old_bbot;
  if ( display_degree <= 2 )
    _p.set_DISPLAY_DEGREE (0);

  bool chk = false;
  int  m   = static_cast<int>(old_bbot.size());
  for ( int i = 0 ; i < m ; ++i ) {
    if ( old_bbot[i] == _EB_ ) {
      p1_bbot[i] = _OBJ_;
      chk = true;
    }
    else if ( old_bbot[i] == _OBJ_ )
      p1_bbot[i] = _NOTHING_;
  }

  if ( !chk )
    return false;

  _p.set_F_TARGET       ( 0.0     );
  _p.set_BB_OUTPUT_TYPE ( p1_bbot );
  _p.reset_stats_file();
  _p.check ( false );

  // modify evaluator:
  Evaluator           * old_ev = _ev_control.get_evaluator();
  Phase_One_Evaluator * p1ev   = new Phase_One_Evaluator ( _p , _out , *old_ev );
  _ev_control.set_evaluator ( p1ev );

  // run MADS with modified parameters:
  stop_reason = run();

  // reset stopping condition:
  if ( stop_reason == _F_TARGET_REACHED_ ) {
    stop        = false;
    stop_reason = _NO_STOP_;
  }
  else
    stop = true;

  // restore evaluator:
  _ev_control.set_evaluator ( old_ev );
  delete p1ev;

  // restore parameters:
  _p.set_F_TARGET       ( old_f_target   );
  _p.set_BB_OUTPUT_TYPE ( old_bbot       );
  _p.set_STATS_FILE     ( old_stats_file_name , old_stats_file );
  _p.set_DISPLAY_DEGREE ( display_degree );
  _p.check ( false );

  // reset mesh index:
  Mesh::reset_mesh_index();
    
  // counters:
  _stats.add_p1_iterations ( _stats.get_iterations() - old_it  );
  _stats.add_p1_bbe        ( _stats.get_bb_eval   () - old_bbe );

  // compute true values of f and h for best feas./inf. points:
  const Eval_Point * x = _barrier.get_best_feasible();
  if ( x ) {
    Eval_Point * modifiable_x = &Cache::get_modifiable_point ( *x );
    old_ev->compute_f ( *modifiable_x );
    old_ev->compute_h ( *modifiable_x );
    _barrier.reset();
    if ( _barrier.insert(*x) != _FULL_SUCCESS_ )
      return false;
  }
  else
    return false;

  if ( display_degree > 2 )
    _out << endl << "end of Mads phase one" << endl;

  return true;
}

/*---------------------------------------------------------*/
/*                    one MADS iteration                   */
/*---------------------------------------------------------*/
void Mads::iteration ( bool              & stop               ,
		       stop_type         & stop_reason        ,
		       Direction         & feas_success_dir   ,
		       Direction         & infeas_success_dir ,
		       const Eval_Point *& new_feas_inc       ,
		       const Eval_Point *& new_infeas_inc       ) {

  success_type success = _UNSUCCESSFUL_;

  // force quit (by pressing ctrl-c):
  if ( !stop && Mads::_force_quit ) {
    stop        = true;
    stop_reason = _CTRL_C_;
    return;
  }

  // poll center selection:
  _barrier.select_poll_center ( success );

  // displays at the beginning of an iteration:
  if ( _p.get_display_degree() > 2 )
    display_iteration_begin ();

  // SEARCH:
  // -------
  search ( stop , stop_reason , success ,
	   new_feas_inc , new_infeas_inc );

  // POLL:
  // -----
  if ( success != _FULL_SUCCESS_ )
    poll ( stop               ,
	   stop_reason        ,
	   success            ,
	   feas_success_dir   ,
	   infeas_success_dir ,
	   new_feas_inc       ,
	   new_infeas_inc       );

  // UPDATES:
  // --------
  if (!stop) {

    // mesh update:
    int old_l = Mesh::get_mesh_index();
    Mesh::update ( success );
    _std_signature.check_min_mesh_sizes ( _p                     ,
					  Mesh::get_mesh_index() ,
					  stop                   ,
					  stop_reason              );

    // display:
    if ( _p.get_updates_display_degree() > 2 )
      _out << endl
	   << "Mesh update:" << endl
	   << "\tprevious mesh index: " << old_l << endl
	   << "\tnew mesh index     : " << Mesh::get_mesh_index()
	   << endl << endl;

    if (stop)
      Mesh::set_mesh_index ( old_l );

    // periodic cache save (overwrite=false):
    if ( _p.get_cache_save_period() > 0 &&
	 _stats.get_iterations()%_p.get_cache_save_period() == _p.get_cache_save_period() - 1 )
      _ev_control.save_caches ( false );
  }

  // number of iterations:
  _stats.add_iteration();
  if ( !stop                       &&
       _p.get_max_iterations() > 0 &&
       _stats.get_iterations() >= _p.get_max_iterations() ) {
    stop        = true;
    stop_reason = _MAX_ITER_REACHED_;
  }

  // displays at the end of an iteration:
  if ( _p.get_display_degree() > 2 ) {
    _out << endl << "end of iteration " << _stats.get_iterations()-1 << endl;
    display_iteration_end ( _p.get_display_degree() ,
			    stop                    ,
			    stop_reason             ,
			    success                 ,
			    new_feas_inc            ,
			    new_infeas_inc            );
  }
}

/*---------------------------------------------------------*/
/*                        the poll                         */
/*---------------------------------------------------------*/
void Mads::poll ( bool              & stop               ,
		  stop_type         & stop_reason        ,
		  success_type      & success            ,
		  Direction         & feas_success_dir   ,
		  Direction         & infeas_success_dir ,
		  const Eval_Point *& new_feas_inc       ,
		  const Eval_Point *& new_infeas_inc       ) {

  if (stop)
    return;

  int display_degree = _p.get_poll_display_degree();
  success            = _UNSUCCESSFUL_;
  new_feas_inc       = NULL;
  new_infeas_inc     = NULL;

  if ( display_degree > 2 )
    _out << endl << "iteration " << _stats.get_iterations()
	 << ", begin of Mads::poll()" << endl;

  set<Priority_Eval_Point>        poll_pts; // the list of poll trial points:
  list<Direction>::const_iterator it , end;
  const Signature               * signature;
  Eval_Point                    * pt;
  const Eval_Point              * x;
  const Direction               * dir;
  int                             i , k;
  poll_center_type                pc_type;
  int                             n;
  int                             m          = _p.get_bb_nb_outputs();
  const vector<bb_input_type>   & bbit       = _p.get_bb_input_types();
  int                             offset     = 0;
  int                             mesh_index = Mesh::get_mesh_index();

  // poll centers:
  const Eval_Point * poll_centers[2] , * poll_center;
  poll_centers[0] = _barrier.get_poll_center();
  poll_centers[1] = (_p.use_sec_poll_center()) ?
    _barrier.get_sec_poll_center() :
    NULL;

  if ( !poll_centers[0] && !poll_centers[1] )
    throw NOMAD::Exception ( __FILE__ , __LINE__ ,
			     "Mads::poll(): could not get a poll center" );
  // success directions:
  const Direction * pfeas_success_dir
    = feas_success_dir.is_defined  () ? &feas_success_dir   : NULL;
  const Direction * pinfeas_success_dir
    = infeas_success_dir.is_defined() ? &infeas_success_dir : NULL;

  // loop on the two poll centers:
  // -----------------------------
  poll_type i_pc = _PRIMARY_;
  poll_center    = poll_centers[_PRIMARY_];

  while ( true ) {

    if ( poll_center ) {

      if ( display_degree > 3 ) {
	_out << endl << "\t";
	if ( i_pc == _SECONDARY_ )
	  _out << "secondary ";
	_out << "poll center: ( ";
	poll_center->Point::display ( _out , " " , 2 , Point::get_display_limit() );
	_out << ")\n" << endl;
      }

      // get the poll center's signature:
      signature = poll_center->get_signature();
      if ( !signature )
	throw NOMAD::Exception ( __FILE__ , __LINE__ ,
	      "Mads::poll(): could not get the poll center's signature" );

      n = signature->get_n();
      if ( n != poll_center->size() )
	throw NOMAD::Exception ( __FILE__ , __LINE__ ,
	      "Mads::poll(): the poll center has an incompatible signature" );
      
      // get directions from the signature:
      list<Direction> dirs;
      signature->get_directions ( dirs , i_pc , Mesh::get_mesh_index() );
      
      if ( !stop && dirs.empty() ) {
	if ( display_degree > 2 )
	  _out << "Mads::poll(): could not get directions: stop" << endl;
	stop        = true;
	stop_reason = _TRIAL_PT_CONSTRUCTION_;
	return;
      }

      // poll center type:
      pc_type = ( poll_center->is_feasible ( _p.get_h_min() ) ) ?
	_FEASIBLE_ : _INFEASIBLE_;

      end = dirs.end();

      // displays:
      if ( display_degree > 3 ) {
	_out << "\n\tpoll directions (include mesh size parameter):\n";
	k = offset;
	for ( it = dirs.begin() ; it != end ; ++it )
	  _out << "\t\tdir " << k++ << " : " << *it << "\n";
      }

      int           nb_dirs = static_cast<int>(dirs.size());
      Random_Pickup rp ( nb_dirs );

      // creation of the poll trial points:
      k = 0;
      for ( it = dirs.begin() ; it != end ; ++it ) {
	
	it->set_index ( offset + k );
	
	dir = &(*it);  
	
	pt = new Eval_Point;
	pt->set                  ( n , m       );
	pt->set_signature        ( signature   );
	pt->set_direction        ( dir         );
	pt->set_mesh_index       ( &mesh_index );
	pt->set_poll_center_type ( pc_type     );
	  
	// random direction ?
	if ( dir_is_random ( dir->get_type() ) )
	  pt->set_rand_eval_priority ( rp.pickup() );

	// pt = poll_center + dir: with a particular case for binary variables
	// equal to 1 with dir=1: the variables are set to 0 (1+1=0 in binary):
	for ( i = 0 ; i < n ; ++i )
	  (*pt)[i] =
	    ( bbit[i] == _BINARY_ && (*dir)[i] == 1.0 && (*poll_center)[i] == 1.0 ) ?
	    0.0 : (*pt)[i] = (*poll_center)[i] + (*dir)[i];

	// we check that the new poll trial point is different than the poll center
	// (this happens when the mesh size becomes near epsilon):
	if ( !stop && pt->Point::operator == ( *poll_center ) ) {
	  delete pt;
	  if ( display_degree > 2 )
	    _out << "Mads::poll(): could not generate poll trial points: stop" << endl;
	  stop        = true;
	  stop_reason = _TRIAL_PT_CONSTRUCTION_;
	  return;
	}
	  
	// add the new point to the ordered list of poll trial points:
	add_eval_point ( poll_pts                ,
			 pt                      ,
			 display_degree          ,
			 _p.get_snap_to_bounds() ,
			 pfeas_success_dir       ,
			 pinfeas_success_dir     ,
			 NOMAD::Double()         ,
			 NOMAD::Double()           );
	
	++k;
      }

      offset = nb_dirs;
	
      // display the re-ordered list of poll trial points:
      if ( display_degree > 3 ) {
	_out << "\n\tre-ordered list of " << static_cast<int>(poll_pts.size())
	     << " poll trial points:\n";
	set<Priority_Eval_Point>::const_iterator end2 = poll_pts.end();
	for ( set<Priority_Eval_Point>::const_iterator it2 = poll_pts.begin() ;
	      it2 != end2 ; ++it2 ) {
	  x =  it2->get_point();
	  _out << "\t\tpoint #" << x->get_tag() << " : (";
	  x->Point::display ( _out , " " , 2 , Point::get_display_limit() );
	  _out << ")";
	  if ( x->get_direction() )
	    _out << " (dir " << x->get_direction()->get_index() << ")";
	  _out << endl;
	}
	_out << endl;
      }
    }
   
    // loop increment:
    if ( i_pc == _PRIMARY_ ) {
      i_pc        = _SECONDARY_;
      poll_center = poll_centers[_SECONDARY_];
    }
    else
      break;
  }

  const Point * old_inf_inc = _barrier.get_best_infeasible();

  _stats.add_poll_pts ( static_cast<int>(poll_pts.size()) );

  // eval_list_of_points (2/4, poll):
  // --------------------------------
  eval_list_of_points ( _POLL_              ,
			poll_pts            ,
			stop                ,
			stop_reason         ,
			new_feas_inc        ,
			new_infeas_inc      ,
			success             ,
			pfeas_success_dir   ,
			pinfeas_success_dir   );

  if ( success == _FULL_SUCCESS_ )
    _stats.add_poll_success();

  // success directions:
  if ( new_feas_inc && new_feas_inc->get_direction() )
    feas_success_dir = *(new_feas_inc->get_direction());
  else
    feas_success_dir.clear();

  if ( new_infeas_inc && new_infeas_inc->get_direction() )
    infeas_success_dir = *(new_infeas_inc->get_direction());
  else
    infeas_success_dir.clear();

  // check the PEB constraints: if we have a new best infeasible incumbent from
  // another infeasible incumbent:
  if ( _p.get_barrier_type() == _PEB_P_ &&
       new_infeas_inc                   &&
       new_infeas_inc->get_direction()  &&
       old_inf_inc                      &&
       *old_inf_inc == *new_infeas_inc - *new_infeas_inc->get_direction() )
                    // *new_infeas_inc - *new_infeas_inc->get_direction() == poll center
    _barrier.check_PEB_constraints ( *new_infeas_inc , display_degree );
  
  // final display:
  if ( display_degree > 2 )
    _out << endl << "iteration " << _stats.get_iterations()
	 << ", end of Mads::poll()" << endl;
}

/*---------------------------------------------------------*/
/*                        the search                       */
/*---------------------------------------------------------*/
void Mads::search ( bool              & stop           ,
		    stop_type         & stop_reason    ,
		    success_type      & success        ,
		    const Eval_Point *& new_feas_inc   ,
		    const Eval_Point *& new_infeas_inc   ) {

  int display_degree = _p.get_search_display_degree();
  success            = _UNSUCCESSFUL_;

  if ( !_p.get_speculative_search() )
    return;

  // first display:
  if ( display_degree > 2 )
    _out << endl << "iteration " << _stats.get_iterations()
	 << ", begin of Mads::search()" << endl;

  // speculative search:
  if ( _p.get_speculative_search() )
    speculative_search ( display_degree ,
			 stop           ,
			 stop_reason    ,
			 success        ,
			 new_feas_inc   ,
			 new_infeas_inc   );

  // Latin-Hypercube (LH) search:
  if ( success != _FULL_SUCCESS_ && _p.get_LH_search_pi() > 0 )
    LH_search ( false          ,
		_std_signature ,
		stop           ,
		stop_reason    ,
		success        ,
		new_feas_inc   ,
		new_infeas_inc   );

  if ( success == _FULL_SUCCESS_ ) {
    if ( display_degree > 2 )
      _out << endl << "iteration " << _stats.get_iterations()
	   << ", end of Mads::search()" << endl;
    return;
  }

  // final display:
  if ( display_degree > 2 )
    _out << endl << "iteration " << _stats.get_iterations()
	 << ", end of Mads::search()" << endl;
}

/*-----------------------------------------------------------*/
/*              MADS Latin-Hypercube (LH) search             */
/*                (based on a given signature)               */
/*-----------------------------------------------------------*/
void Mads::LH_search ( bool                initial_search ,
		       const Signature   & signature      ,
		       bool              & stop           ,
		       stop_type         & stop_reason    ,
		       success_type      & success        ,
		       const Eval_Point *& new_feas_inc   ,
		       const Eval_Point *& new_infeas_inc   ) {
  if (stop)
    return;

  int display_degree = _p.get_search_display_degree();

  if ( display_degree > 3 )
    _out << "\n\t" << ( initial_search ? "initial " : "")
	 << "Latin-Hypercube search\n\n";
  
  int   i;
  int   p = initial_search ? _p.get_LH_search_p0() : _p.get_LH_search_pi();
  int   n = signature.get_n();

  set<Priority_Eval_Point> search_pts;
  Eval_Point             * x;
  int                      m          = _p.get_bb_nb_outputs();
  int                      mesh_index = Mesh::get_mesh_index();
  int                      pm1        = p-1;

  // mesh sizes:
  Point delta_m_max (n);
  signature.get_delta_m ( delta_m_max , Mesh::get_min_mesh_index() );

  Point delta_m;
  if ( !initial_search )
    signature.get_delta_m ( delta_m , mesh_index );

  const Point & fixed_variables = signature.get_fixed_variables();
  
  // pts contains n points of dimension p: each of these points contains
  // p different values for each variable:
  Point ** pts = new Point * [n];

  // creation of p search points:
  for ( int k = 0 ; k < p ; ++k ) {

    x = new Eval_Point;
    x->set ( n , m );
    x->set_signature  ( &signature  );
    x->set_mesh_index ( &mesh_index );

    for ( i = 0 ; i < n ; ++i ) {

      if ( k==0 ) {
	if ( fixed_variables[i].is_defined() )
	  pts[i] = new Point ( p , fixed_variables[i] );
	else {
	  pts[i] = new Point(p);
	  LH_values_for_var_i ( i , p , !initial_search , delta_m , delta_m_max , *pts[i] );
	}
      }

      (*x)[i] = (*pts[i])[k];

      if ( k == pm1 )
	delete pts[i];
    }

    if ( display_degree > 3 ) {
      _out << "\tLH search point #" << x->get_tag()
	   << ": (";
      x->Point::display ( _out , " " , 2 , Point::get_display_limit() );
      _out << ")" << endl;
    }

    // add the new point to the ordered list of search trial points:
    add_eval_point ( search_pts      ,
		     x               ,
		     display_degree  ,
		     false           ,  // snap_to_bounds = false
		     NULL            ,
		     NULL            ,
		     NOMAD::Double() ,
		     NOMAD::Double()   );
  }

  delete [] pts;

  _stats.add_LH_pts ( static_cast<int>(search_pts.size()) );

  // eval_list_of_points (4/4, LH search):
  // -------------------------------------
  eval_list_of_points ( _LH_SEARCH_  ,
			search_pts     ,
			stop           ,
			stop_reason    ,
			new_feas_inc   ,
			new_infeas_inc ,
			success        ,
			NULL           ,
			NULL             );

  if ( success == _FULL_SUCCESS_ )
    _stats.add_LH_success();
}

/*-----------------------------------------------------------*/
/*  LH search: decide p values for one variable (index ind)  */
/*-----------------------------------------------------------*/
/*     (if no bounds, values are scaled with the largest     */
/*      delta_m value obtained so far)                       */
/*-----------------------------------------------------------*/
void Mads::LH_values_for_var_i ( int                   ind         ,
				 int                   p           ,
				 bool                  project     ,
				 const Point         & delta_m     ,
				 const Point         & delta_m_max ,
				 Point               & x             ) const {

  bb_input_type         bbit  = _p.get_bb_input_types()[ind];
  const NOMAD::Double & lb    = _p.get_lb()            [ind];
  const NOMAD::Double & ub    = _p.get_ub()            [ind];
  const NOMAD::Double & dmax  = delta_m_max            [ind];

  NOMAD::Double         delta = project ? delta_m[ind] : NOMAD::Double();

  // caterogical variables: not yet coded:
  if ( bbit == _CATEGORICAL_ )
    return;

  Random_Pickup rp(p);
  int           i;
  bool          rounding  = (bbit != _CONTINUOUS_);
  bool          lb_def    = lb.is_defined();
  bool          ub_def    = ub.is_defined();
  double        w         = ( ( lb_def && ub_def ) ? ub.value()-lb.value() : 1.0 ) / p;
  NOMAD::Double v;

  // main loop:
  for ( i = 0 ; i < p ; ++i ) {

    // both bounds exist:
    if ( lb_def && ub_def )
      v = lb + ( i + rand()/D_INT_MAX ) * w;

    // one of the bounds does not exist:
    else {

      // lb exists, and ub not: mapping [0;1] --> [lb;+INF[
      if ( lb_def )
	v = lb + 10 *
	  dmax * sqrt ( - log ( DEFAULT_EPSILON + ( i + rand()/D_INT_MAX ) * w ) );

      // lb does not exist:
      else {

	// ub exists, and lb not: mapping [0;1] --> ]-INF;ub]
	if ( ub_def )
	  v = ub - dmax * 10 *
	    sqrt ( -log ( DEFAULT_EPSILON + ( i + rand()/D_INT_MAX ) * w ) );
	
	// there are no bounds: mapping [0;1] --> ]-INF;+INF[
	else
	  v = (rand()%2 ? -1.0 : 1.0) * dmax * 10 *
	    sqrt ( - log ( DEFAULT_EPSILON + ( i + rand()/D_INT_MAX ) * w ) );
      }
    }

    // rounding:
    if ( rounding )
      v = v.round();

    // projection to mesh:
    if ( project ) {
      v = (v / delta).round() * delta;
      if ( v > ub )
	v = ub;
      else if ( v < lb )
	v = lb;
    }

    // affectation + permutation:
    x[rp.pickup()] = v;
  }
}

/*-----------------------------------------------------------*/
/*                   MADS speculative search                 */
/*-----------------------------------------------------------*/
/*  x_k = x_{k-1} + \Delta^m_{k-1} d                         */
/*  s_k = x_{k-1} +        \Delta^m_{k-1} d  if l_{k-1} > 0  */
/*            or    \tau * \Delta^m_{k-1} d  else            */
/*-----------------------------------------------------------*/
/*  the directions that we use already contain \Delta^m:     */
/*     dir = \Delta^m_{k-1} d                                */
/*-----------------------------------------------------------*/
void Mads::speculative_search ( int                 display_degree ,
				bool              & stop           ,
				stop_type         & stop_reason    ,
				success_type      & success        ,
				const Eval_Point *& new_feas_inc   ,
				const Eval_Point *& new_infeas_inc   ) {

  if ( stop || (!new_feas_inc && !new_infeas_inc) )
    return;

  if ( display_degree > 3 )
    _out << "\n\tspeculative search" << endl;

  set<Priority_Eval_Point> search_pts;
  int                      lkm1;  // l_{k-1}
  int                      lk;    // l_k
  int                      n;
  const Signature        * signature;
  Point                    delta_m_k;
  Point                    delta_m_km1;
  Point                    factor;
  Point                    xkm1;
  Eval_Point             * sk;
  const Eval_Point       * x[2];
  x[0] = new_feas_inc;
  x[1] = new_infeas_inc;
  
  for ( int i = 0 ; i < 2 ; ++i ) {

    if ( x[i]                             &&
	 x[i]->get_direction()            &&
	 x[i]->get_direction()->is_mads() &&
	 x[i]->get_mesh_index()              ) {

      // get the x_k's signature:
      signature = x[i]->get_signature();
      if ( !signature )
	throw NOMAD::Exception ( __FILE__ , __LINE__ ,
	      "Mads::speculative_search(): could not get the signature" );

      xkm1  = *x[i] - *x[i]->get_direction();

      lk = lkm1 = *x[i]->get_mesh_index();
      Mesh::update ( _FULL_SUCCESS_ , lk );

      n = signature->get_n();
      delta_m_k   = n;
      delta_m_km1 = n;

      signature->get_delta_m ( delta_m_k   , lk   );
      signature->get_delta_m ( delta_m_km1 , lkm1 );

      // multiplicative factor: takes into account the fact that
      // the direction contains \Delta^m_k:
      try {
	factor = delta_m_k / delta_m_km1;
      }
      catch ( NOMAD::Double::Invalid_Value & ) {
	if ( display_degree > 2 )
	  _out << "could not compute speculative search point: stop" << endl;
	stop        = true;
	stop_reason = _TRIAL_PT_CONSTRUCTION_;
	return;
      }

      if ( lkm1 <= 0 )
	factor *= Mesh::get_mesh_update_basis();
   
      // speculative search point:
      Direction new_dir ( n , 0.0 , x[i]->get_direction()->get_type() );
      new_dir.Point::operator = ( factor * *x[i]->get_direction() );

      sk = new Eval_Point;
      sk->set ( n , _p.get_bb_nb_outputs() );
      sk->set_signature  ( signature );
      sk->set_direction  ( &new_dir );
      sk->set_mesh_index ( &lk );

      sk->Point::operator = ( xkm1 + new_dir );

      if ( display_degree > 3 ) {
	_out << "\tspeculative search point #" << sk->get_tag()
	     << ": (";
	sk->Point::display ( _out ," " , 2 , Point::get_display_limit() );
	_out << ")" << endl;
      }

      // add the new point to the ordered list of search trial points:
      add_eval_point ( search_pts              ,
		       sk                      ,
		       display_degree          ,
		       _p.get_snap_to_bounds() ,
		       NULL                    ,
		       NULL                    ,
		       NOMAD::Double()         ,
		       NOMAD::Double()           );
    }
  }

  _stats.add_spec_pts ( static_cast<int>(search_pts.size()) );

  // eval_list_of_points (3/4, speculative search):
  // ----------------------------------------------
  eval_list_of_points ( _MADS_SEARCH_  ,
			search_pts     ,
			stop           ,
			stop_reason    ,
			new_feas_inc   ,
			new_infeas_inc ,
			success        ,
			NULL           ,
			NULL             );

  if ( success == _FULL_SUCCESS_ )
    _stats.add_spec_success();
}

/*-----------------------------------------------------------*/
/*  add an Eval_Point to the list of points to be evaluated  */
/*          (the point is also snapped to bounds)            */
/*-----------------------------------------------------------*/
void Mads::add_eval_point ( set<Priority_Eval_Point>   & pts                ,
			    Eval_Point                *& x                  ,
			    int                          display_degree     ,
			    bool                         snap_to_bounds     ,
			    const Direction           *  feas_success_dir   ,
			    const Direction           *  infeas_success_dir ,
			    const NOMAD::Double       &  f_sgte             ,
			    const NOMAD::Double       &  h_sgte               ) const {
 
  // snap to bounds:
  if ( snap_to_bounds && x->snap_to_bounds() ) {
    if ( display_degree > 3 ) {
      _out << "\n\tpoint #" << x->get_tag() << " ";
      if ( x->get_direction() && x->get_direction()->get_index() >= 0 )
	_out << "(dir " << x->get_direction()->get_index() << ") ";
      _out << "has been snapped to bounds" << endl;
    }

    if ( x->get_direction() && x->get_direction()->norm() == 0.0 ) {
      
      if (display_degree > 3 )
	_out << "\tpoint #" << x->get_tag() << " is flushed (||dir||==0)" << endl;
      delete x;
      x = NULL;

      return;
    }
  }

  // creation of the Priority_Eval_Point:
  Priority_Eval_Point pep ( x , _p.get_h_min() );

  // ordering elements of Priority_Eval_Point's:
  // -------------------------------------------

  // 1. user criterion: TODO

  // 2. surrogate values for f and h:
  pep.set_f_sgte ( f_sgte );
  pep.set_h_sgte ( h_sgte );
  
  // 3. angle with last successful directions:
  if ( x->get_direction() ) {
    
    if ( feas_success_dir && x->get_poll_center_type() == _FEASIBLE_ )
      pep.set_angle_success_dir ( feas_success_dir->get_angle ( *x->get_direction() ) );

    if ( infeas_success_dir && x->get_poll_center_type() == _INFEASIBLE_ )
      pep.set_angle_success_dir ( infeas_success_dir->get_angle ( *x->get_direction() ) );
  }

  // 4. angle with simplex gradient:
  // TODO

  // insertion of the point in pts:
  // ------------------------------
  size_t size_before = pts.size();

  // insertion:
  pts.insert(pep);

  // do not delete this point at the end of this method:
  if ( pts.size() > size_before )
    pep.do_not_delete_point();
}

/*---------------------------------------------------------*/
/*                          x0 eval                        */
/*---------------------------------------------------------*/
void Mads::eval_x0 (  bool      & stop        ,
		      stop_type & stop_reason   ) {
  const vector<Point *> & x0s           = _p.get_x0s();
  const string          & x0_cache_file = _p.get_x0_cache_file();
  if ( x0s.empty() && x0_cache_file.empty() )
    return;

  if ( _p.get_display_degree() > 2 )
    _out << endl << "begin of Mads::eval_x0()" << endl;

  Eval_Point                             * pt;
  set<Priority_Eval_Point>                 pts;
  set<Priority_Eval_Point>::const_iterator it , end;
  size_t                                   k;
  int                                      m   = _p.get_bb_nb_outputs();
  int                                      n   = _p.get_dimension();
  ostringstream                            err;

  // x0s from vector Parameters::_x0s:
  // ---------------------------------
  size_t x0s_size = x0s.size();
  for ( k = 0 ; k < x0s_size ; ++k ) {
    
    // the current starting point has to be in dimension n:
    if ( x0s[k]->size() != n ) {
      err << "starting point " << *x0s[k] << " is not of dimension " << n;
      throw NOMAD::Exception ( __FILE__ , __LINE__ , err.str() );
    }

    // creation of the Eval_Point:
    pt = new Eval_Point;
    pt->set           ( *x0s[k] , m     );
    pt->set_signature ( &_std_signature );

    add_eval_point    ( pts                     ,
			pt                      ,
			_p.get_display_degree() ,
			false                   ,
			NULL                    ,
			NULL                    ,
			NOMAD::Double()         ,
			NOMAD::Double()           );
  }

  // x0 from a cache file:
  // ---------------------
  if ( !x0_cache_file.empty() ) {

    Cache            & cache = _ev_control.get_cache();
    const Eval_Point * x;

    // another cache file (this file won't be modified):
    if ( x0_cache_file != _p.get_cache_file() ) {

      Cache  x0_cache ( _out , _p.get_display_degree() , false );

      string file_name = _p.get_problem_dir() + x0_cache_file;

      if ( !x0_cache.load ( file_name ) ) {
 	err << "could not load (or create) the cache file " << file_name;
 	throw NOMAD::Exception ( __FILE__ , __LINE__ , err.str() );
      }

      // we copy all the temporary cache points into pts:
      x = x0_cache.begin();
      while (x) {

	pt = new Eval_Point;
	pt->set( *x , m );

	if ( x->get_signature() )
	  pt->set_signature ( x->get_signature() );
	else if ( x->size() == n )
	  pt->set_signature ( &_std_signature );

	if ( pt->get_signature() )
	  add_eval_point ( pts                     ,
			   pt                      ,
			   _p.get_display_degree() ,
			   false                   ,
			   NULL                    ,
			   NULL                    ,
			   NOMAD::Double()         ,
			   NOMAD::Double()           );
	else {
	  if ( _p.get_display_degree() > 1 )
	    _out << "\nWarning (" << __FILE__ << ", " << __LINE__
		 << "): could not add the evaluation point " << *pt
		 << " (no signature)\n\n";
	  delete pt;
	}
	x = x0_cache.next();
      }

      // insertion of this temporary cache in the algorithm's cache:
      cache.insert ( x0_cache );
    }

    // x0 cache file and the algorithm's cache file are the same:
    else {
      x = cache.begin();
      while (x) {
	pt = &cache.get_modifiable_point ( *x );

	if ( x->get_signature() )
	  pt->set_signature ( x->get_signature() );
	else if ( x->size() == n )
	  pt->set_signature ( &_std_signature );

	if ( pt->get_signature() )
	  add_eval_point ( pts                     ,
			   pt                      ,
			   _p.get_display_degree() ,
			   false                   ,
			   NULL                    ,
			   NULL                    ,
			   NOMAD::Double()         ,
			   NOMAD::Double()           );
	else {
	  if ( _p.get_display_degree() > 1 )
	    _out << "\nWarning (" << __FILE__ << ", " << __LINE__
		 << "): could not add the evaluation point " << *pt
		 << " (no signature)\n\n";
	  delete pt;
	}
	x = cache.next();
      }
    }
  }

  // display of all starting points:
  if ( _p.get_display_degree() > 3 ) {
    k = 0;
    _out << "\nx0 eval points:\n";
    end = pts.end();
    for ( it = pts.begin() ; it != end ; ++it ) {
      _out << "\tpoint " << static_cast<int>(k++) << ": ( ";
      it->get_point()->Point::display ( _out , " " , 2 , Point::get_display_limit() );
      _out << " )\n";
    }
  }

  success_type       success;   
  const Eval_Point * new_feas_inc   = NULL;
  const Eval_Point * new_infeas_inc = NULL;

  // eval_list_of_points (1/4, x0 eval.):
  // ------------------------------------
  eval_list_of_points ( _X0_           ,
			pts            ,
			stop           ,
			stop_reason    ,
			new_feas_inc   ,
			new_infeas_inc ,
			success        ,
			NULL           ,
			NULL             );
  
  if ( !stop &&
       ( success==_UNSUCCESSFUL_ ||
	 ( _p.get_barrier_type() == _EB_ && !_barrier.get_best_feasible() ) ) ) {
    stop        = true;
    stop_reason = _X0_EVAL_;
  }

  // displays:
  display_iteration_end ( _p.get_display_degree() ,
			  stop                    ,
			  stop_reason             ,
			  success                 ,
			  new_feas_inc            ,
			  new_infeas_inc            );

  // stop the algorithm if no iterations are allowed:
  if ( !stop && _p.get_max_iterations() == 0 ) {
    stop        = true;
    stop_reason = _MAX_ITER_REACHED_;
  }

  if ( _p.get_display_degree() > 2 ) 
    _out << endl << "end of Mads::eval_x0()" << endl;
}

/*-----------------------------------------------------------------*/
/*  use the Evaluator_Control object to evaluate a list of points  */
/*-----------------------------------------------------------------*/
void Mads::eval_list_of_points ( search_type                search             ,
				 set<Priority_Eval_Point> & points             ,
				 bool                     & stop               ,
				 stop_type                & stop_reason        ,
				 const Eval_Point        *& new_feas_inc       ,
				 const Eval_Point        *& new_infeas_inc     ,
				 success_type             & success            ,
				 const Direction         *  feas_success_dir   ,
				 const Direction         *  infeas_success_dir   ) {

  list<const Eval_Point *> evaluated_pts;

  // use the surrogate to sort the eval. points ?
  if ( _p.get_sgte_eval_sort() && points.size() > 1 ) {

    // define all points as surrogates:
    for ( set<Priority_Eval_Point>::iterator it = points.begin() ;
	  it != points.end() ; ++it )
      Cache::get_modifiable_point ( *it->get_point() ).set_surrogate ( true );

    // evaluate the surrogates:
    _ev_control.eval_list_of_points ( search         ,
				      points         ,
				      NULL           , // barrier = NULL
				      stop           ,
				      stop_reason    ,
				      new_feas_inc   ,
				      new_infeas_inc ,
				      success        ,
				      evaluated_pts    );
    if (stop)
      return;

    // construct a new list of trial points that will be
    // ordered using surrogate values:
    Eval_Point * x;
    int          display_degree = (search==_POLL_) ? _p.get_poll_display_degree  () :
                                                     _p.get_search_display_degree() ;

    list<const Eval_Point *>::const_iterator end = evaluated_pts.end();
    for ( list<const Eval_Point *>::const_iterator it2 = evaluated_pts.begin() ;
	  it2 != end ; ++it2 ) {

      // the current point is still a surrogate (normal case): we create
      // a new Eval_Point point:
      if ( (*it2)->is_surrogate() ) {
	
	// Eval_Point construction:
	x = new Eval_Point;
	x->set ( (*it2)->size() , _p.get_bb_nb_outputs() );
	x->set_signature  ( (*it2)->get_signature () );
	x->set_direction  ( (*it2)->get_direction () );
	x->set_mesh_index ( (*it2)->get_mesh_index() );
	x->Point::operator = ( **it2 );

	// display:
	if ( display_degree > 3 ) {
	  _out << "\n\tnew " << search << " point #" << x->get_tag()
	       << ": (";
	  x->Point::display ( _out , " " , 2 , Point::get_display_limit() );
	  _out << ")" << endl;
	}

	// add the new point to the ordered list of  trial points:
	add_eval_point ( points                  ,
			 x                       ,
			 display_degree          ,
			 _p.get_snap_to_bounds() ,
			 feas_success_dir        ,
			 infeas_success_dir      ,
			 (*it2)->get_f()         ,
			 (*it2)->get_h()           );
      }

      // the current point is not anymore a surrogate: it is directly added
      // in the new list of eval_points
      // (it corresponds to the particular case where a surrogate point
      //  was found in the 'true' cache, see Evaluator_Control::eval_lop()):
      else {

	x = &Cache::get_modifiable_point(**it2);

	// display:
	if ( display_degree > 3 ) {
	  _out << "\n\t" << search << " point #" << x->get_tag()
	       << ": (";
	  x->Point::display ( _out , " " , 2 , Point::get_display_limit() );
	  _out << ")" << endl;
	}

	add_eval_point ( points             ,
			 x                  ,
			 0                  ,  // display_degree = 0
			 false              ,  // snap_to_bounds = false
			 feas_success_dir   ,
			 infeas_success_dir ,
			 NOMAD::Double()    ,
			 NOMAD::Double()      );
      }

    }
  }

  if (stop)
    return;

  // evaluate the list of points on the 'true' function:
  _ev_control.eval_list_of_points ( search         ,
				    points         ,
				    &_barrier      ,
				    stop           ,
				    stop_reason    ,
				    new_feas_inc   ,
				    new_infeas_inc ,
				    success        ,
				    evaluated_pts    );
}

/*---------------------------------------------------------*/
/*                         MADS display                    */
/*---------------------------------------------------------*/
void Mads::display ( ostream & out ) const {

  int display_degree = _p.get_display_degree();
  if ( display_degree <= 1 )
    return;

  // detailed display:
  if ( _p.get_display_degree() > 2 ) {
    _out << endl
	 << "\tcache:\n\t------\n" << _ev_control.get_cache() << endl
	 << "\tconstraints handling:\n\t---------------------\n" << _barrier << endl
	 << "\tstats:\n\t------\n";
    _stats.display ( _out , _p.get_display_degree() );
    _out << endl
	 << "\tmisc:\n\t-----" << endl
	 << "\tmesh index            : min="
	 << Mesh::get_min_mesh_index() << ", max="
	 << Mesh::get_max_mesh_index() << ", last="
	 << Mesh::get_mesh_index() << endl
	 << "\tbest feasible point   : ";
    if ( _barrier.get_best_feasible() ) {
      _out << "( ";
      _barrier.get_best_feasible()->Point::display ( _out , " " , -1 , -1 );
      _out << " ) h=" << _barrier.get_best_feasible()->get_h()
	   << " f=" << _barrier.get_best_feasible()->get_f() << endl;
    }
    else
      _out << "NULL\n";
    _out << "\tbest infeasible point : ";
    if ( _barrier.get_best_infeasible() ) {
      _out << "( ";
      _barrier.get_best_infeasible()->Point::display ( _out , " " , -1 , -1 );
      _out << " ) h=" << _barrier.get_best_infeasible()->get_h()
	   << " f=" << _barrier.get_best_infeasible()->get_f() << endl;
    }
    else
      _out << "NULL\n";
    _out << endl;
  }

  // normal display:
  else {

    // black-box evaluations:
    _out << endl
	 << "black-box evaluations : " << _stats.get_bb_eval() << endl;

    // output stats:
    if ( _stats.get_stat_sum().is_defined() )
      _out << "stat sum              : " << _stats.get_stat_sum() << endl;
    if ( _stats.get_stat_avg().is_defined() )
      _out << "stat avg              : " << _stats.get_stat_avg() << endl;

    // best solutions:
    if ( !_barrier.get_best_feasible() && !_barrier.get_best_infeasible() )
      _out << "no solution" << endl;
    else {

      if ( _barrier.get_best_infeasible() ) {
	_out << "best infeasible point : (";
	_barrier.get_best_infeasible()->Point::display ( _out , " " , -1 , -1 );
	_out << ") h=" << _barrier.get_best_infeasible()->get_h()
	     << " f="  << _barrier.get_best_infeasible()->get_f() << endl;
      }

      if ( _barrier.get_best_feasible() ) {
	_out << "best feasible point   : (";
	_barrier.get_best_feasible()->Point::display ( _out , " " , -1 , -1 );
	_out << ") h=" << _barrier.get_best_feasible()->get_h()
	     << " f="  << _barrier.get_best_feasible()->get_f() << endl;
      }
    }
  }
}

/*---------------------------------------------------------*/
/*  display mesh and poll sizes for a given signature and  */
/*  the current mesh index                                 */
/*---------------------------------------------------------*/
void Mads::display_deltas ( const Signature & s ) const {
  Point delta_m , delta_p;
  s.get_delta_m ( delta_m , Mesh::get_mesh_index() );
  s.get_delta_p ( delta_p , Mesh::get_mesh_index() );
  _out << "\tmesh size             : ( " << delta_m << " )\n"
       << "\tpoll size             : ( " << delta_p << " )"
       << endl;
}

/*---------------------------------------------------------*/
/*         displays at the beginning of an iteration       */
/*---------------------------------------------------------*/
void Mads::display_iteration_begin ( void ) const {
  const Signature * s1 = NULL;
  _out << endl
       << "iteration " << _stats.get_iterations() << ":\n"
       << "--------------\n"
       << "\tblack-box evaluations : " << _stats.get_bb_eval() << endl
       << "\tbest feasible point   : ";
  if ( _barrier.get_best_feasible() ) {
    _out << "( ";
    _barrier.get_best_feasible()->Point::display ( _out , " " , 2 , Point::get_display_limit() );
    _out << " ) h=" << _barrier.get_best_feasible()->get_h()
	 << " f=" << _barrier.get_best_feasible()->get_f()
	 << endl;
  }
  else
    _out << "NULL\n";
  _out << "\tbest infeasible point : ";
  if ( _barrier.get_best_infeasible() ) {
    _out << "( ";
    _barrier.get_best_infeasible()->Point::display ( _out , " " , 2 , Point::get_display_limit() );
    _out << " ) h=" << _barrier.get_best_infeasible()->get_h()
	 << " f=" << _barrier.get_best_infeasible()->get_f()
	 << endl;
  }
  else
    _out << "NULL\n";

  _out << "\tpoll center           : ";
  if ( _barrier.get_poll_center() ) {
    _out << "( ";
    _barrier.get_poll_center()->Point::display ( _out , " " , 2 , Point::get_display_limit() );
    _out << " )\n";

    s1 = _barrier.get_poll_center()->get_signature();
    if (s1)
      display_deltas ( *s1 );
  }
  else
    _out << "NULL\n";

  if ( _barrier.get_sec_poll_center() ) {
    _out << "\tsec. poll center      : ( ";
    _barrier.get_sec_poll_center()->Point::display ( _out , " " , 2 , Point::get_display_limit() );
    _out << " )\n";
    const Signature * s2 = _barrier.get_sec_poll_center()->get_signature();
    if ( s2 && (!s1 || s1 != s2) )
      display_deltas ( *s2 );
  }
  _out << "\tmesh index            : " << Mesh::get_mesh_index() << endl
       << "\th_max                 : " << _barrier.get_h_max()   << endl;
}

/*---------------------------------------------------------*/
/*             displays at the end of an iteration         */
/*---------------------------------------------------------*/
void Mads::display_iteration_end ( int                display_degree ,
				   bool               stop           ,
				   stop_type          stop_reason    ,
				   success_type       success        ,
				   const Eval_Point * new_feas_inc   ,
				   const Eval_Point * new_infeas_inc   ) const {

  if ( display_degree > 2 ) {
    _out << endl
	 << "\tterminate MADS           : ";
    display_yes_or_no ( stop , _out );
    _out << endl;
    if (stop)
      _out << "\ttermination cause        : " << stop_reason << endl;
    _out << "\titeration status         : " << success << endl;
    _out << "\tnew feasible incumbent   : ";
    if (new_feas_inc)
      _out << *new_feas_inc << endl;
    else
      _out << "NULL\n";
    _out << "\tnew infeasible incumbent : ";
    if (new_infeas_inc)
      _out << *new_infeas_inc << endl;
    else
      _out << "NULL\n";
  }
}
