/*-------------------------------------------------------------------------------------*/
/*  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 "Directions.hpp"
using namespace std;

/*---------------------------------------------------------*/
/*                         destructor                      */
/*---------------------------------------------------------*/
Directions::~Directions ( void ) {
  if ( _lt_initialized )
    for ( int i = 0 ; i <= 2*L_LIMITS ; ++i )
      delete _bl[i];
  delete [] _primes;
}

/*---------------------------------------------------------*/
/*                      other initializations              */
/*---------------------------------------------------------*/
void Directions::init ( const set<direction_type> & dir_types          ,
			const set<direction_type> & sec_poll_dir_types ,
			bb_output_type              barrier_type         ) {

  // default values for _dir_types:
  if ( _dir_types.empty() )
    _dir_types = dir_types;

  // default values for _sec_poll_dir_types:
  if ( barrier_type == _PB_ || barrier_type == _PEB_P_ ) {
    if ( _sec_poll_dir_types.empty() )
      _sec_poll_dir_types = sec_poll_dir_types;
    else if ( _sec_poll_dir_types.find(_NO_DIRECTION_) != _sec_poll_dir_types.end() )
      _sec_poll_dir_types.clear();
  }
  else
    _sec_poll_dir_types.clear();

  // is_orthomads: true if at least one direction is of type Ortho-MADS:
  _is_orthomads = false;
  set<direction_type>::const_iterator it , end = _dir_types.end();
  for ( it = _dir_types.begin() ; it != end ; ++it )
    if ( dir_is_orthomads (*it) ) {
      _is_orthomads = true;
      break;
    }
  if ( !_is_orthomads ) {
    end = _sec_poll_dir_types.end();
    for ( it = _sec_poll_dir_types.begin() ; it != end ; ++it )
      if ( dir_is_orthomads (*it) ) {
	_is_orthomads = true;
	break;
      }
  }
}

/*---------------------------------------------------------*/
/*                  LT-MADS initializations                */
/*---------------------------------------------------------*/
void Directions::lt_mads_init ( void ) {
  for ( int i = 0 ; i <= 2*L_LIMITS ; ++i ) {
    _bl   [i] = NULL;
    _hat_i[i] = -1;
  }
  _lt_initialized = true;
}

/*---------------------------------------------------------*/
/*                   Ortho-MADS initializations            */
/*---------------------------------------------------------*/
void Directions::ortho_mads_init ( int halton_seed ) {

  _is_orthomads = true;

  if (!_primes) {
    _primes = new int[_nc];
    construct_primes( _nc , _primes );
  }

  _halton_seed = ( halton_seed < 0 ) ? _primes[_nc-1] : halton_seed;

#ifdef DEBUG
  _out << endl << "Ortho-MADS Halton seed (t0) = " << _halton_seed << endl << endl;
#endif

  _max_halton_index = _halton_seed;
}

/*---------------------------------------------------------*/
/*                         set_binary                      */
/*---------------------------------------------------------*/
void Directions::set_binary ( void ) {

  _binary       = true;
  _is_orthomads = false;

  _dir_types.clear();
  _dir_types.insert(_GPS_BINARY_);

  if ( !_sec_poll_dir_types.empty() ) {
    _sec_poll_dir_types.clear();
    _sec_poll_dir_types.insert(_GPS_BINARY_);
  }
}

/*---------------------------------------------------------*/
/*                set_categorical (not yet coded)          */
/*---------------------------------------------------------*/
void Directions::set_categorical ( void ) {
  _categorical  = true;
  _is_orthomads = false;
  _dir_types.clear();
  _sec_poll_dir_types.clear();
}

/*----------------------------------------------------------------------*/
/*  compute binary directions (only if all group-variables are binary)  */
/*----------------------------------------------------------------------*/
void Directions::compute_binary_directions ( list<Direction> & dirs ,
					     poll_type         poll   ) const {
  // _GPS_BINARY_ n directions:
  Direction * pd;
  for ( int i = 0 ; i < _nc ; ++i ) {
    dirs.push_back ( Direction ( _nc , 0.0 , _GPS_BINARY_ ) );
    pd = &(*(--dirs.end()));
    (*pd)[i] = 1.0;
  }
}

/*---------------------------------------------------------*/
/*            get the directions for a given mesh          */
/*---------------------------------------------------------*/
void Directions::compute ( list<Direction> & dirs         ,
			   poll_type         poll         ,
			   int               mesh_index   ,
			   int               halton_index   ) {
  dirs.clear();

  // categorical variables: not yet coded: no directions are computed:
  if ( _categorical )
    return;

  // binary variables: we use special directions:
  if ( _binary ) {
    compute_binary_directions ( dirs , poll );
    return;
  }

  NOMAD::Double     pow_gps , cst;
  const Direction * bl;
  Direction       * pd;
  int               i , j , k , hat_i;
  list<Direction>::const_iterator it2 , end2;

  /*-----------------------------------*/
  /*  loop on the types of directions  */
  /*-----------------------------------*/
  const set<direction_type> & dir_types = (poll==_PRIMARY_) ?
    _dir_types : _sec_poll_dir_types;

  set<direction_type>::const_iterator end = dir_types.end() ;
  for ( set<direction_type>::const_iterator it = dir_types.begin() ;
	it != end ; ++it ) {

    if ( *it == _UNDEFINED_DIRECTION_ || *it == _NO_DIRECTION_ )
      continue;

    /*--------------*/
    /*  Ortho-MADS  */
    /*--------------*/
    if ( dir_is_orthomads (*it) ) {

      // Ortho-MADS initializations:
      if (!_primes)
	ortho_mads_init ( _halton_seed );

      // halton index:
      if ( halton_index < 0 ) {
	
	halton_index = ( mesh_index >= _max_mesh_index ) ?
	  mesh_index + _halton_seed :
	  _max_halton_index + 1;
  
	if ( halton_index > _max_halton_index )
	  _max_halton_index = halton_index;
      }
      
      Direction     halton_dir ( _nc , 0.0 , *it );
      NOMAD::Double alpha_t_l;

      if ( compute_halton_dir ( halton_index ,
				mesh_index   ,
				alpha_t_l    ,
				halton_dir     ) ) {

#ifdef DEBUG
	_out << "\tcompute Ortho-MADS directions with:\n"
	     << "\t\tHalton index (tk) = " << halton_index << endl
	     << "\t\tmesh index   (lk) = " << mesh_index   << endl
	     << "\t\talpha     (tk,lk) = " << alpha_t_l    << endl
	     << "\t\tHalton direction  = (";
	halton_dir.Point::display ( _out , " " , -1 , -1 );
	_out << " )" << endl;
#endif

	// Ortho-MADS 2n:
	// --------------
	if ( *it == _ORTHO_2N_ ) {

	  // creation of the 2n directions:
	  Direction ** H = new Direction * [2*_nc];
	  for ( i = 0 ; i < _nc ; ++i ) {
	    dirs.push_back ( Direction ( _nc , 0.0 , _ORTHO_2N_ ) );
	    H[i    ] = &(*(--dirs.end()));
	    dirs.push_back ( Direction ( _nc , 0.0 , _ORTHO_2N_ ) );
	    H[i+_nc] = &(*(--dirs.end()));
	  }

	  // Householder transformations on the 2n directions:
	  householder ( halton_dir , H );

	  delete [] H;
	}

	// Ortho-MADS 1 or Ortho-MADS 2:
	// -----------------------------
	else {
	  dirs.push_back ( halton_dir );
	  if ( *it == _ORTHO_2_ )
	    dirs.push_back ( -halton_dir );
	}
      }
    }

    /*-----------*/
    /*  LT-MADS  */
    /*-----------*/
    else if ( dir_is_ltmads (*it) ) {

      if ( !_lt_initialized)
	lt_mads_init();

      bl = get_bl ( mesh_index , *it , hat_i );
      
      // LT-MADS 1 or LT-MADS 2: -b(l) and/or b(l):
      // ------------------------------------------
      if ( *it == _LT_1_ || *it == _LT_2_ ) {

	dirs.push_back ( - *bl );
	if ( *it == _LT_2_ )
	  dirs.push_back ( *bl );
	
      }

      // LT-MADS 2n or LT-MADS n+1:
      // --------------------------
      else {

	// row permutation vector:
	int * row_permutation_vector  = new int [_nc];
	row_permutation_vector[hat_i] = hat_i;

	Random_Pickup rp(_nc);

	for ( i = 0 ; i < _nc ; ++i )
	  if ( i != hat_i ) {
	    j = rp.pickup();
	    if ( j == hat_i )
	      j = rp.pickup();
	    row_permutation_vector[i] = j;
	  }
 
	rp.reset();

	for ( j = 0 ; j < _nc ; ++j ) {

	  i = rp.pickup();

	  if ( i != hat_i ) {

	    dirs.push_back ( Direction ( _nc , 0.0 , *it ) );
	    pd = &(*(--dirs.end()));

	    create_lt_direction ( mesh_index , *it , i , hat_i , pd );

	    permute_coords ( *pd , row_permutation_vector );
	  }
	  else
	    dirs.push_back ( *bl );

	  // completion to maximal basis:
	  if ( *it == _LT_2N_ )
	    dirs.push_back ( Direction ( - *(--dirs.end()) , _LT_2N_ ) );
	}
	
	delete [] row_permutation_vector;

	// completion to minimal basis:
	if ( *it == _LT_NP1_ ) {

	  dirs.push_back ( Direction ( _nc , 0.0 , _LT_NP1_ ) );
	  pd = &(*(--dirs.end()));

	  end2 = --dirs.end();
	  for ( it2 = dirs.begin() ; it2 != end2 ; ++it2 ) {
	    for ( i = 0 ; i < _nc ; ++i )
	      (*pd)[i] -= (*it2)[i];
	  }
	}
      }
    }

    /*-------*/
    /*  GPS  */
    /*-------*/
    else {

      // GPS for binary variables: should'nt be here:
      if ( *it == _GPS_BINARY_ ) {
#ifdef DEBUG
	_out << "\nWarning (" << __FILE__ << ", " << __LINE__
	     << "): tried to generate binary directions at wrong the place)\n\n";
#endif
	dirs.clear();
	return;
      }

      // this value represents the non-zero values in GPS directions
      // (it is tau^{|lk|/2}, and it is such that Delta_k^m * pow_gps = Delta_k^p):
      if ( !pow_gps.is_defined() )
	pow_gps = pow ( Mesh::get_mesh_update_basis() , abs(mesh_index) / 2.0 );

      // GPS 2n, static:
      // ---------------
      if ( *it == _GPS_2N_STATIC_ ) {

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

	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_2N_STATIC_ ) );
	  pd = &(*(--dirs.end()));
	  (*pd)[i] = pow_gps;

	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_2N_STATIC_ ) );
	  pd = &(*(--dirs.end()));
	  (*pd)[i] = -pow_gps;
	}
      }

      // GPS 2n, random:
      // ---------------
      else if ( *it == _GPS_2N_RAND_ ) {

	int v , cnt;

	list  <int>::const_iterator end3;
	list  <int>::iterator       it3;
	list  <int> rem_cols;
	vector<int> rem_col_by_row(_nc);

	// creation of the 2n directions:
	vector<Direction *> pdirs (2*_nc);

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

	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_2N_RAND_ ) );
	  pdirs[i] = &(*(--dirs.end()));
	  (*pdirs[i])[i] = pow_gps;

	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_2N_RAND_ ) );
	  pdirs[i+_nc] = &(*(--dirs.end()));
	  (*pdirs[i+_nc])[i] = -pow_gps;

	  rem_cols.push_back(i);
	  rem_col_by_row[i] = i;
	}
      
	// random perturbations:
	for ( i = 1 ; i < _nc ; ++i ) {
     	  if ( rem_col_by_row[i] > 0 ) { 
	    v = rand()%3 - 1; // v in { -1 , 0 , 1 }
	    if ( v ) {

	      // we decide a (i,j) couple:
	      k    = rand()%(rem_col_by_row[i]);
	      cnt  = 0;
	      end3 = rem_cols.end();
	      it3  = rem_cols.begin();
	      while ( cnt != k ) {
		++it3;
		++cnt;
	      }
	      j = *it3;

	      // the perturbation:
	      (*pdirs[i])[j] = (*pdirs[j+_nc])[i] = -v * pow_gps;
	      (*pdirs[j])[i] = (*pdirs[i+_nc])[j] =  v * pow_gps;
	      
	      // updates:
	      rem_cols.erase(it3);      
	      it3 = rem_cols.begin();
	      while ( *it3 != i )
		++it3;
	      rem_cols.erase(it3);
	      for ( k = i+1 ; k < _nc ; ++k )
		rem_col_by_row[k] -= j<k ? 2 : 1;
	    }       
	  }
	}

      }
      
      // GPS n+1, static:
      // ----------------
      else if ( *it == _GPS_NP1_STATIC_ ) {

	// (n+1)^st direction:
	dirs.push_back ( Direction ( _nc , -pow_gps , _GPS_NP1_STATIC_ ) );

	// first n directions:
	for ( i = 0 ; i < _nc ; ++i ) {
	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_STATIC_ ) );
	  pd = &(*(--dirs.end()));
	  (*pd)[i] = pow_gps;
	}
      }

      // GPS n+1, random:
      // ----------------
      else if ( *it == _GPS_NP1_RAND_ ) {

	// (n+1)^st direction:
	dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_RAND_ ) );
	Direction * pd1 = &(*(--dirs.end()));
	
	NOMAD::Double d;

	// first n directions:
	for ( i = 0 ; i < _nc ; ++i ) {
	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_RAND_ ) );
	  pd = &(*(--dirs.end()));
	  
	  d = rand()%2 ? -pow_gps : pow_gps;
	  (*pd )[i] =  d;
	  (*pd1)[i] = -d;
	}
      }

      // GPS n+1, static, uniform angles:
      // --------------------------------
      else if ( *it == _GPS_NP1_STATIC_UNIFORM_ ) {

	cst = pow_gps * sqrt(static_cast<double>(_nc)*(_nc+1))/_nc;

	// n first directions:
	for ( j = _nc-1 ; j >= 0 ; --j ) {
	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_STATIC_UNIFORM_ ) );
	  pd = &(*(--dirs.end()));
	  
	  k = _nc-j-1;

	  for ( i = 0 ; i < k ; ++i )
	    (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));

	  (*pd)[k] = (cst * (j+1)) / sqrt(static_cast<double>(j+1)*(j+2));

	}

	// (n+1)^st direction:
	dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_STATIC_UNIFORM_ ) );
	pd = &(*(--dirs.end()));
	for ( i = 0 ; i < _nc ; ++i )
	  (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));

      }

      // GPS n+1, random, uniform angles:
      // --------------------------------
      // (we apply the procedure defined in "Pattern Search Methods for user-provided points:
      // application to molecular geometry problems", by Alberto, Nogueira, Rocha and Vicente,
      // SIOPT 14-4, 1216-1236, 2004, doi:10.1137/S1052623400377955)
      else if ( *it == _GPS_NP1_RAND_UNIFORM_ ) {

	cst = pow_gps * sqrt(static_cast<double>(_nc)*(_nc+1))/_nc;

	// n first directions:
	for ( j = _nc-1 ; j >= 0 ; --j ) {
	  dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_STATIC_UNIFORM_ ) );
	  pd = &(*(--dirs.end()));

	  k = _nc-j-1;

	  for ( i = 0 ; i < k ; ++i )
	    (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));

	  (*pd)[k] = (cst * (j+1)) / sqrt(static_cast<double>(j+1)*(j+2));

	}

	// (n+1)^st direction:
	dirs.push_back ( Direction ( _nc , 0.0 , _GPS_NP1_STATIC_UNIFORM_ ) );
	pd = &(*(--dirs.end()));
	for ( i = 0 ; i < _nc ; ++i )
	  (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));

	// random perturbations:
	Random_Pickup             rp   (_nc);
	vector<bool>              done (_nc);
	bool                      chg_sgn;
	list<Direction>::iterator it;
	NOMAD::Double             tmp;

	end2 = dirs.end();
	for ( i = 0 ; i < _nc ; ++i )
	  done[i] = false;

	for ( i = 0 ; i < _nc ; ++i ) {
    
	  k = rp.pickup();

	  if ( i != k && !done[i] ) {

	    chg_sgn = ( rand()%2 != 0 );

	    for ( it = dirs.begin() ; it != end2 ; ++it ) {
	      tmp      = (*it)[i];
	      (*it)[i] = ( chg_sgn ? -1.0 : 1.0 ) * (*it)[k];
	      (*it)[k] = tmp;
	    }

	    done[i] = done[k] = true;
	  }
	  else
	    if ( rand()%2 )
	      for ( it = dirs.begin() ; it != end2 ; ++it )
		(*it)[i] = -(*it)[i];

	}

      }

    }

  } // end of loop on the types of directions

  // _max_mesh_index update:
  if ( mesh_index > _max_mesh_index )
    _max_mesh_index = mesh_index;
}

/*---------------------------------------------------------*/
/*            compute the halton direction q(t,l)          */
/*---------------------------------------------------------*/
bool Directions::compute_halton_dir ( int             halton_index ,
				      int             mesh_index   ,
				      NOMAD::Double & alpha_t_l    ,
				      Direction     & halton_dir     ) const {

  alpha_t_l.clear();

  int           i;
  NOMAD::Double d , norm = 0.0;
  Point         b (_nc);

  for ( i = 0 ; i < _nc ; ++i ) {
    d = get_phi ( halton_index , _primes[i] );
    b[i] = 2*d - 1.0;
    norm += b[i].pow2();
  }
  norm = norm.sqrt();

  // desired squared norm for q:
  NOMAD::Double target = ( halton_dir.get_type() == _ORTHO_2N_ ) ?
    pow ( Mesh::get_mesh_update_basis() , abs(mesh_index) / 2.0 ) :
    pow ( Mesh::get_mesh_update_basis() , abs(mesh_index)       );

  NOMAD::Double x  = target.sqrt();
  NOMAD::Double fx = eval_ortho_norm ( x , norm , target , b , halton_dir );
  NOMAD::Double y  = 1.0 , fy , delta_x , abs_dx , min , delta_min , diff , eps;
  bool          inc , possible , set_eps = false;
  int           cnt = 0;

  while ( fx != target ) {

    // safety test:
    if ( ++cnt > 1000 ) {
#ifdef DEBUG
      _out << "\nWarning (" << __FILE__ << ", " << __LINE__
	   << "): could not compute Halton direction for (t="
	   << halton_index << ",l=" << mesh_index << ")\n\n";
#endif
      return false;
    }

    if ( set_eps ) {
      eps     = 1e-8;
      set_eps = false;
    }
    else
      eps = 0.0;

    inc = ( fx < target );
     
    possible = false;
    min      = 1e+20;
    for ( i = 0 ; i < _nc ; ++i ) {
      
      if ( b[i] != 0.0 ) {
	
	if ( b[i] > 0.0 ) {
	  if ( inc )
	    diff  =  0.5+eps;
	  else
	    diff  = -0.5-eps;
	  
	}
	else {
	  if ( inc )
	    diff  = -0.5-eps;
	  else
	    diff  =  0.5+eps;
	}
	
	delta_x = norm * ( halton_dir[i] + diff) / b[i] - x;
	
	y = x + delta_x;

	if ( y > 0 ) {
	  abs_dx = delta_x.abs();
	  if ( abs_dx < min ) {
	    min       = abs_dx;
	    delta_min = delta_x;
	    possible  = true;
	  }
	}
      }
    }

    if ( !possible ) {
      string err =
	"Directions::compute_halton_dir(): could not compute a Halton direction (n="
	+ NOMAD::itos (_nc         )
	+ " t=" + NOMAD::itos (halton_index)
	+ " l=" + NOMAD::itos (mesh_index  ) + ")";
      throw NOMAD::Exception ( __FILE__ , __LINE__ , err );
    }

    y  = x + delta_min;
    fy = eval_ortho_norm ( y , norm , target , b , halton_dir );

    if ( fx == fy ) {
      set_eps = true;
      continue;
    }
    
    if ( fy==target )
      break;
    
    if ( inc && fy > target && fx > 0 ) {
      eval_ortho_norm ( x , norm , target , b , halton_dir );
      break;
    }
    
    if ( !inc && fy < target && fy > 0 )
      break;
    
    x  = y;
    fx = fy;

  }

  alpha_t_l = y;

  return true;
}

/*---------------------------------------------------------------*/
/*  compute the squared norm of normalized(2u_t-e) (Ortho-MADS)  */
/*---------------------------------------------------------------*/
NOMAD::Double Directions::eval_ortho_norm ( const NOMAD::Double & x      ,
					    const NOMAD::Double & norm   ,
					    const NOMAD::Double & target ,
					    const Point         & b      ,
					    Point               & new_b    ) const {
  NOMAD::Double fx = 0.0;

  for ( int i = 0 ; i < _nc ; ++i ) {
    new_b[i] = ( x * b[i] / norm ).round();
    fx += new_b[i]*new_b[i];
  }
  
  return fx;
}

/*--------------------------------------------------------*/
/*  get the expression of an integer t in inverse base p  */
/*--------------------------------------------------------*/
NOMAD::Double Directions::get_phi ( int t , int p ) const {
  int   div;
  int  size = int ( ceil ( log(static_cast<double>(t+1)) / log(static_cast<double>(p)) ) );
  int    ll = t;
  NOMAD::Double d = 0.0;
  for ( int i = 0 ; i < size ; ++i ) {
    div = int ( pow ( static_cast<double>(p) , size-i-1 ) );
    d += ( ll / div ) * pow ( static_cast<double>(p) , i-size );
    ll = ll % div;
  }
  return d;
}

/*----------------------------------------------------------------*/
/*  . Householder transformation to generate _nc directions the   */
/*    Halton direction                                            */
/*  . compute also H[i+nc] = -H[i] (completion to 2n directions)  */
/*----------------------------------------------------------------*/
void Directions::householder ( const Direction & halton_dir , Direction ** H ) const {
  int i , j;

  NOMAD::Double norm2 = halton_dir.squared_norm() , v , h2i;

  for ( i = 0 ; i < _nc ; ++i ) {
    
    h2i = 2 * halton_dir[i];

    for ( j = 0 ; j < _nc ; ++j ) {

      // H[i]:
      (*H[i    ])[j] =  v = (i==j) ? norm2 - h2i * halton_dir[j] : - h2i * halton_dir[j];

      // -H[i]:
      (*H[i+_nc])[j] = -v;
    }
  }
}

/*---------------------------------------------------------*/
/*               get the LT-MADS b(l) direction            */
/*---------------------------------------------------------*/
const Direction * Directions::get_bl ( int            mesh_index ,
				       direction_type dtype      ,
				       int          & hat_i        ) {

  Direction * bl = _bl    [ mesh_index + L_LIMITS ];
  hat_i          = _hat_i [ mesh_index + L_LIMITS ];

  if ( !bl ) {
    hat_i = -1;
    create_lt_direction ( mesh_index , dtype , -1 , hat_i , bl );
  }

  return bl;
}

/*---------------------------------------------------------*/
/*              create a new LT-MADS direction             */
/*---------------------------------------------------------*/
/*           (if hat_i == -1, a new b(l) direction         */
/*            is created and hat_i is set)                 */
/*---------------------------------------------------------*/
void Directions::create_lt_direction ( int                    mesh_index ,
				       direction_type         dtype      ,
				       int                    diag_i     ,
				       int                  & hat_i      ,
				       Direction           *& dir          ) {
  int i_pow_tau =
    static_cast<int>
    ( ceil ( pow ( Mesh::get_mesh_update_basis() , abs(mesh_index) / 2.0 ) ) );

  int j = diag_i+1;

  // new b(l) direction:
  if ( hat_i < 0 ) {
    _hat_i [ mesh_index + L_LIMITS ] = diag_i = hat_i = rand()%_nc;
    _bl    [ mesh_index + L_LIMITS ] = dir = new Direction ( _nc , 0.0 , dtype );

    j = 0;
  }

  (*dir)[diag_i] = (rand()%2) ? -i_pow_tau : i_pow_tau;

  for ( int k = j ; k < _nc ; ++k )
    if ( k != hat_i ) {
      (*dir)[k] = rand()%i_pow_tau;
      if ( rand()%2 && (*dir)[k] > 0.0 )
	(*dir)[k] = -(*dir)[k];
    }

#ifdef DEBUG
  if ( j==0 )
    _out << "new LT-MADS b(l) direction: b(" << mesh_index << ") = "
	 << *dir << endl << endl;
#endif
}

/*---------------------------------------------------------*/
/*            permute the coordinates of a direction       */
/*---------------------------------------------------------*/
void Directions::permute_coords ( Direction & dir , int * permutation_vector ) const {
  Point tmp = dir;
  for ( int i = 0 ; i < _nc ; ++i )
    dir [ permutation_vector[i] ] = tmp[i];
}

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

  out << "n=" << _nc << ", types={ ";
  set<direction_type>::const_iterator it , end = _dir_types.end();
  for ( it = _dir_types.begin() ; it != end ; ++it )
    out << "[" << *it << "] ";
  out << "}";
  out << ", sec_poll_types={ ";
  end = _sec_poll_dir_types.end();
  for ( it = _sec_poll_dir_types.begin() ; it != end ; ++it )
    out << "[" << *it << "] ";
  out << "}";

  if ( _is_orthomads ) {
    out << ", halton_seed=";
    if ( _halton_seed > 0 )
      out << _halton_seed;
    else
      out << "auto";
  }
}
