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

/*-----------------------------------*/
/*   static members initialization   */
/*-----------------------------------*/
#ifdef DEBUG
int Point::_cardinality     = 0;
int Point::_max_cardinality = 0;
#endif

int Point::_display_limit   = DEFAULT_POINT_DISPLAY_LIMIT;

/*-----------------------------------------------------------*/
/*                         constructor                       */
/*-----------------------------------------------------------*/
Point::Point ( int n , const NOMAD::Double & d ) : _n (n) , _coords (NULL) {
#ifdef DEBUG
  ++Point::_cardinality;
  if ( Point::_cardinality > Point::_max_cardinality )
    ++Point::_max_cardinality;
#endif
  if (_n > 0) {
    _coords = new NOMAD::Double [_n];  
    if ( d.is_defined() )
      fill ( _coords , _coords+_n , d );
  }
  else
    _n = 0;
}

/*-----------------------------------------------------------*/
/*                        copy constructor                   */
/*-----------------------------------------------------------*/
Point::Point ( const Point & p ) : _n (p._n) , _coords (NULL) {
#ifdef DEBUG
  ++Point::_cardinality;
  if ( Point::_cardinality >= Point::_max_cardinality )
    ++Point::_max_cardinality;
#endif
  if ( _n > 0 ) {
    NOMAD::Double       * p1 =   _coords = new NOMAD::Double [_n];
    const NOMAD::Double * p2 = p._coords;
    for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 )
      *p1 = *p2;
  }
}

/*-----------------------------------------------*/
/*                    destructor                 */
/*-----------------------------------------------*/
Point::~Point (void) {
#ifdef DEBUG
  clear();
  --Point::_cardinality;
#endif
  delete [] _coords;
}

/*-----------------------------------------------*/
/*   This method changes the point's dimension   */
/*   and sets all values to d                    */
/*-----------------------------------------------*/
void Point::reset ( int n , const NOMAD::Double & d )  {
  if ( n <= 0 ) {
    _n = 0; 
    delete [] _coords;
    _coords = NULL;
  }
  else {
    if ( _n != n ) {
      delete [] _coords;
      _n      = n;
      _coords = new NOMAD::Double [_n];
    }
    if ( d.is_defined() )
      fill ( _coords , _coords+_n , d );
  }
}

/*----------------------------------------------------------------*/
/*  This method changes the array's dimension (and keeps values)  */
/*----------------------------------------------------------------*/
void Point::resize ( int n ) {
  if ( n == _n )
    return;
  if ( n <= 0 ) {
    _n = 0; 
    delete [] _coords;
    _coords = NULL;
    return;
  }
  NOMAD::Double * new_coords = new NOMAD::Double [n];
  if ( _coords ) {

    int min = ( n < _n ) ? n : _n;
    
    NOMAD::Double       * p1 = new_coords;
    const NOMAD::Double * p2 = _coords;

    for ( int i = 0 ; i < min ; ++i , ++p1 , ++p2 )
      *p1 = *p2;

    delete [] _coords;
  }
  _coords = new_coords;
  _n      = n;
}

/*-----------------------------------------------------------*/
/*                       '[]' operators                      */
/*-----------------------------------------------------------*/

// const version:
const NOMAD::Double & Point::operator [] ( int i ) const {
  if ( !_coords )
    throw Not_Defined ( __FILE__ , __LINE__ ,
			"operator x[i] (const): 'x' not defined" );
  if ( i < 0 || i >= _n )
    throw Bad_Access ( __FILE__ , __LINE__ ,
		       "operator x[i] (const): 'i' outside the array's bounds." );
  return _coords[i];
}

// non-const version:
NOMAD::Double & Point::operator [] ( int i ) {
  if ( !_coords )
    throw Not_Defined ( __FILE__ , __LINE__ ,
			"operator x[i]: 'x' not defined" );
  if ( i < 0 || i >= _n )
    throw Bad_Access ( __FILE__ , __LINE__ ,
		       "operator x[i] (const): 'i' outside the array's bounds." );
  return _coords[i];
}

/*-----------------------------------------------------------*/
/*                     affectation operator                  */
/*-----------------------------------------------------------*/
const Point & Point::operator = ( const Point & p ) {
  if ( this == &p )
    return *this;
	
  if ( _n != p._n ) {
    delete [] _coords;
    _n = p._n;
    if (_n > 0)
      _coords = new NOMAD::Double [_n];
    else
      _coords = NULL;
  }

  NOMAD::Double       * p1 =   _coords;
  const NOMAD::Double * p2 = p._coords;
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 )
    *p1 = *p2;
	
  return *this;
}

/*-----------------------------------------------------------*/
/*                             display                       */
/*-----------------------------------------------------------*/
void Point::display ( ostream & out , const string & sep , int w , int lim ) const {

  int nm1 = _n-1;
  
  // for a limited display of maximum lim elements:
  if ( lim > 0 && lim < _n ) {

    int l1 = (lim + 1) / 2 , l2 = lim / 2 , i;

    // first coordinates:
    for ( i = 0 ; i < l1 ; ++i )
      out << setw(w) << _coords[i] << sep;

    // separator:
    out << "..." << sep;

    // last coordinates:
    for ( i = _n - l2 ; i < nm1 ; ++i )
      out << setw(w) << _coords[i] << sep;
  }

  // normal display (lim <= 0 or lim >= _n):
  else {
    const NOMAD::Double * p = _coords;
    for ( int i = 0 ; i < nm1 ; ++i , ++p )
      out << setw(w) << *p << sep;
  }

  // last coordinate (different because there is no separator after that):
  if ( _n > 0 )
    out << setw(w) << _coords[nm1];
}

/*-----------------------------------------------------------*/
/*                             input                         */
/*-----------------------------------------------------------*/
istream & operator >> ( istream & in , Point & p ) {
  int n = p.size();
  NOMAD::Double * pt = p._coords;
  for ( int k = 0 ; k < n ; ++k , ++pt )
    in >> *pt;
  if ( in.fail() )
    throw Point::Bad_Input ( __FILE__ , __LINE__ , "in >> x: bad input" );
  return in;
}

/*-----------------------------------------------------------*/
/*  set the point's coordinate with the array 'a' of size n  */
/*  the Point's dimension is changed to n                    */
/*-----------------------------------------------------------*/
void Point::set ( int n , const NOMAD::Double * a ) {
  if ( n <= 0 || !a )
    return;

  if ( _n != n ) {
    delete [] _coords;
    _n      = n;
    _coords = new NOMAD::Double [_n];
  }

  NOMAD::Double * p = _coords;
  for ( int k = 0 ; k < _n ; ++k , ++p , ++a )
    *p = *a;
}

/*-----------------------------------------------------------*/
/*        Check if all the point values are defined          */
/*-----------------------------------------------------------*/
bool Point::is_complete ( void ) const {
  const NOMAD::Double * p = _coords;
  for ( int i = 0 ; i < _n ; ++i , ++p )
    if ( !p->is_defined() )
      return false;
  return true;
}

/*---------------------------------------------------------------*/
/*  Check if at least one value is defined in the _coords array  */
/*---------------------------------------------------------------*/
bool Point::is_defined ( void ) const {
  const NOMAD::Double * p = _coords;
  for ( int i = 0 ; i < _n ; ++i , ++p )
    if ( p->is_defined() )
      return true;
  return false;
}

/*-----------------------------------------------------------*/
/*                           negation                        */
/*-----------------------------------------------------------*/
const Point Point::operator - ( void ) const {
  Point tmp (_n);
  NOMAD::Double       * p1 = tmp._coords;
  const NOMAD::Double * p2 = _coords;
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 )
    *p1 = - *p2;
  return tmp;
}

/*----------------------------------------------------------*/
/*                    scalar multiplication                 */
/*----------------------------------------------------------*/
const Point & Point::operator *= ( const NOMAD::Double & d ) {
  NOMAD::Double * p = _coords;
  for ( int k = 0 ; k < _n ; ++k , ++p )
    *p *= d;
  return *this;
}

/*----------------------------------------------------------*/
/*               multiplication of two points               */
/*----------------------------------------------------------*/
const Point Point::operator * ( const Point & p ) const {
  if ( p.size() != _n )
    throw Point::Bad_Operation ( __FILE__ , __LINE__ ,
				 "x * y: x.size != y.size" );
  Point tmp (_n);
  NOMAD::Double       * p1 = tmp._coords;
  const NOMAD::Double * p2 =     _coords;
  const NOMAD::Double * p3 =   p._coords;
    
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 , ++p3 )
    *p1 = *p2 * *p3;

  return tmp;
}

/*----------------------------------------------------------*/
/*                     division of two points               */
/*----------------------------------------------------------*/
const Point Point::operator / ( const Point & p ) const {
  if ( p.size() != _n )
    throw Point::Bad_Operation ( __FILE__ , __LINE__ ,
				 "x / y: x.size != y.size" );
  Point tmp (_n);
  NOMAD::Double       * p1 = tmp._coords;
  const NOMAD::Double * p2 =     _coords;
  const NOMAD::Double * p3 =   p._coords;
    
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 , ++p3 )
    *p1 = *p2 / *p3;

  return tmp;
}

/*----------------------------------------------------------*/
/*                   addition of two points                 */
/*----------------------------------------------------------*/
const Point Point::operator + ( const Point & p ) const {
  if ( p.size() != _n )
    throw Point::Bad_Operation ( __FILE__ , __LINE__ ,
				 "x + y: x.size != y.size" );
  Point tmp (_n);
  NOMAD::Double       * p1 = tmp._coords;
  const NOMAD::Double * p2 =     _coords;
  const NOMAD::Double * p3 =   p._coords;
    
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 , ++p3 )
    *p1 = *p2 + *p3;

  return tmp;
}

/*----------------------------------------------------------*/
/*                  subtraction of two points               */
/*----------------------------------------------------------*/
const Point Point::operator - ( const Point & p ) const {
  if ( p.size() != _n )
    throw Point::Bad_Operation ( __FILE__ , __LINE__ ,
				 "x - y: x.size != y.size" );
  Point tmp (_n);
  NOMAD::Double       * p1 = tmp._coords;
  const NOMAD::Double * p2 =     _coords;
  const NOMAD::Double * p3 =   p._coords;
    
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 , ++p3 )
    *p1 = *p2 - *p3;

  return tmp;
}

/*--------------------------------------------------------------------------*/
/* Comparison operator '<': it is used to find and store the points         */
/*                           in the cache                                   */
/*--------------------------------------------------------------------------*/
bool Point::operator < ( const Point & p ) const {
  if ( this == &p )
    return false;

  if ( _n < p._n )
    return true;
  if ( _n > p._n )
    return false;

  const NOMAD::Double * p1 =   _coords;
  const NOMAD::Double * p2 = p._coords;

  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 ) {

    if ( *p1 < *p2 )
      return true;

    if ( *p1 > *p2 )
      return false;
  }
  return false;
}

/*-----------------------------------------------------------*/
/*                           operator ==                     */
/*-----------------------------------------------------------*/
bool Point::operator == ( const Point & p ) const {
  if ( this == &p )
    return true;
  if ( p.size() != _n )
    return false;

  const NOMAD::Double * p1 =   _coords;
  const NOMAD::Double * p2 = p._coords;
  for ( int k = 0 ; k < _n ; ++k , ++p1 , ++p2 )
    if ( *p1 != *p2 )
      return false;

  return true;
}

/*-----------------------------------------------------------*/
/*        computation of the angle with another point        */
/*-----------------------------------------------------------*/
const NOMAD::Double Point::get_angle ( const Point & y ) const {

  if ( _n != y.size() )
    return NOMAD::Double();

  NOMAD::Double inner_product = 0.0 , norm_x = 0.0 , norm_y = 0.0;

  const NOMAD::Double * p1 =   _coords;
  const NOMAD::Double * p2 = y._coords;

  for ( int i = 0 ; i < _n ; ++i , ++p1 , ++p2 ) {
    norm_x        += *p1 * *p1;
    norm_y        += *p2 * *p2;
    inner_product += *p1 * *p2;
  }
  
  if ( norm_x == 0.0 || norm_y == 0.0 )
    return NOMAD::Double();

  return acos ( ( inner_product / ( norm_x.sqrt() * norm_y.sqrt() ) ).value() );
}
