/******************************************************************************
* NOMAD -  Nonlinear Optimization for Mixed vAriables and Derivatives.        *
* Copyright (C) 2001-2004  Gilles Couture - Ecole Polytechnique,              *
*                          Charles Audet - Ecole Polytechnique,               *
*		           John Dennis - Rice University.                     *
*                                                                             *
* Author: Gilles Couture                                                      *
* École Polytechnique de Montréal - GERAD                                     *
* C.P. 6079, Succ. Centre-ville, Montréal (Québec) H3C 3A7 Canada             *
* Tel:  (514) 340-6053 x6928                                                  *
* Email: Gilles.Couture@gerad.ca                                              *
* Fax: (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. You can find the              *
* GNU General Public License at:                                              *
*             http://www.gnu.org/copyleft/gpl.html                            *
*                                                                             *
* 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, write to the Free Software                 *
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. *
*                                                                             *
* You can find information on the NOMAD software at:                          *
*             http://www.gerad.ca/NOMAD/                                      *
******************************************************************************/

// constraints.cpp -- Constraints class methods.

#include <cmath>
#include <fstream>
#include <iostream>
using namespace std;
#include "constraints.h"
#include "defines.h"
#include "description.h"
#include "parameters.h"

// Private methods:
// ****************************************************************************

// This method calculates the tolerances within which the constraints are
// considered active. It creates the epsilon and proximity arrays, and
// initialize them.
void Constraints::buildArrays()
{
  // The mesh size, mesh size termination criterion, mesh factor and problem
  // dimension are retrieved.
  double delta = parameters->getPollSize();
  double deltaT = parameters->getPollSizeTerm();
  double factor = parameters->getPollBasis();
  int max_dimension = description->getDimension();

  double bound_diff, diff;
  // The arrays are built.
  epsilon_array = new double[max_dimension];
  proximity_array = new int[max_dimension];

  // The epsilon values for each dimension are calculated in this loop.
  for (int i = 0; i < max_dimension; i++)
    {
      // The first part works even when the bounds are not finite.
      // If the mesh size is a termination criterion...
      if (deltaT > 0.0)
	// we multiply it by 4 (arbitrary value).
	epsilon_array[i] = 4 * deltaT;
      else   // mesh size is not a termination criterion.
	// epsilon_array[i] =  (starting) mesh size * (mesh factor ^ -5).
	epsilon_array[i] = delta * pow(factor, -5.0);

      // The second part works for finite bounds.
      // We check if the bounds are finite.
      if ((lower_bounds[i] != -INFINIT) && (upper_bounds[i] != INFINIT) )
	{
	  // We calculate a value relative to the lower and upper bounds, and
	  // put it in 'diff'.
	  bound_diff = upper_bounds[i] - lower_bounds[i];
	  diff = bound_diff * 0.01;
	  // If this value is greater than the one calculated in the first
	  // part, we keep it.
	  if (diff > epsilon_array[i])
	    epsilon_array[i] = diff;
	}

      // The proximity array is initialized to '0' values.
      proximity_array[i] = 0;
    }

  // The epsilon array is shown.
  if (display_factor >= 5)
    showEpsilonArray();
}

// This method outputs the epsilon array.
void Constraints::showEpsilonArray() const
{
  if (epsilon_array != NULL)
    {
      int max_dimension = description->getDimension();
      cout << "Epsilon array: (";
      for (int i = 0; i < (max_dimension - 1); i++)
	cout << epsilon_array[i] << ", ";
      cout << epsilon_array[max_dimension - 1] << ")\n";
    }
}

// This method outputs the proximity array.
void Constraints::showProximityArray() const
{
  if (proximity_array != NULL)
    {
      int max_dimension = description->getDimension();
      cout << "Proximity array: (";
      for (int i = 0; i < (max_dimension - 1); i++)
	cout << proximity_array[i] << ", ";
      cout << proximity_array[max_dimension - 1] << ")\n";
    }
}

// Public methods:
// ****************************************************************************

// Default constructor. This object is created by Evaluator when the problem
// has bound constraints.
Constraints::Constraints(const Description * const des,
			 const Parameters * const par)
  : description(des), parameters(par)
{
  // These arrays are null at the start.
  epsilon_array = NULL;
  proximity_array = NULL;
  display_factor = 0;

  int max_dimension = description->getDimension();
  // The constraint arrays are created.
  lower_bounds = new double[max_dimension];
  upper_bounds = new double[max_dimension];

  // We need pointers to these arrays to scale the bounds.
  double * low = description->getLowerBounds();
  double * up = description->getUpperBounds();
  double * scales_add = description->getScalesAdd();
  double * scales_mult = description->getScalesMult();

  // The constraint arrays are initialized.
  for (int i = 0; i < max_dimension; i++)
    {
      // 'lower_bounds' and 'upper_bounds' are scaled.
      lower_bounds[i] = (low[i] / scales_mult[i]) + scales_add[i];
      upper_bounds[i] = (up[i] / scales_mult[i]) + scales_add[i];
    }

  // This variable will hold the user's choice of poll directions.
  int choice = parameters->getPollDirs();
  // If the user chose GPS n+1 directions or GPS uniform directions, the
  // epsilon and proximity arrays are built.
  if ((choice == 1) || (choice == 2))
    buildArrays();
}

// Destructor. The arrays are destroyed.
Constraints::~Constraints()
{
  if (epsilon_array != NULL)
    delete [] epsilon_array;
  delete [] lower_bounds;
  delete [] upper_bounds;
  if (proximity_array != NULL)
    delete [] proximity_array;
}

// This method checks if the point at the current iteration is outside or
// inside the bounds. It returns true if it's inside, false otherwise.
bool Constraints::check(const double * coord) const
{
  int max_dimension = description->getDimension();
  // Each coordinate has to be inside the bounds.
  for (int i = 0; i < max_dimension; i++)
    {
      // If the coordinate is outside the bounds...
      if ((coord[i] < lower_bounds[i]) || (coord[i] > upper_bounds[i]) )
	// the method will return 'false'.
	return false;
    }
  return true;
}

// The user wants to continue the current run.
void Constraints::continueRun()
{
  // This variable will hold the user's choice of poll directions.
  int choice = parameters->getPollDirs();

  // If the arrays already exist.
  if ((epsilon_array != NULL) && (proximity_array != NULL))
    {
      // If the user didn't chose GPS n+1 directions or GPS uniform directions,
      // the arrays are freed.
      if (!((choice == 1) || (choice == 2)))
				{
					delete [] epsilon_array;
					delete [] proximity_array;
				}
    }
  else   // The arrays aren't created.
    {
      // If the user chose GPS n+1 directions or GPS uniform directions, the
      // arrays are created.
      if ((choice == 1) || (choice == 2))
				buildArrays();
    }
}

// The values of the proximity vector are updated.
void Constraints::updateProximityArray(const double * coord)
{
  // If the proximity array is null we do nothing; it means that the poll
  // directions don't need to be updated.
  if (proximity_array != NULL)
    {
      int i, max_dimension;
      max_dimension = description->getDimension();
      if (display_factor >= 5)
	{
	  cout << "\nUpdating proximity array...\n";
	  cout << "Received point: ";
	  for (i = 0; i < max_dimension - 1; i++)
	    cout << coord[i] << ", ";
	  cout << coord[max_dimension - 1] << "\n";
	}
      // The proximity array is updated in this loop..
      for (i = 0; i < max_dimension; i++)
	{
	  // If the coordinate is within epsilon of the bounds...
	  if ((coord[i] - epsilon_array[i] < lower_bounds[i]) || 
	      (coord[i] + epsilon_array[i] > upper_bounds[i]) )
	    // ... we give it a '1' value.
	    proximity_array[i] = 1;
	  else   // The coordinate is farther inside than epsilon.
	    // The value is '0'.
	    proximity_array[i] = 0;
	}
      // The proximity array is shown.
      if (display_factor >= 5)
	showProximityArray();
    }
}

// Non-member functions.
// ****************************************************************************

// This method outputs this object's attributes.
ostream & operator<<(ostream & os, const Constraints & cons)
{
  int i, max_dimension;
  max_dimension = cons.description->getDimension();

  // The lower bounds array is shown.
  if (cons.lower_bounds != NULL)
    {
      os << "Scaled lower bounds: (";
      for (i = 0; i < (max_dimension - 1); i++)
	os << cons.lower_bounds[i] << ", ";
      os << cons.lower_bounds[max_dimension - 1] << ")\n";
    }
  // The upper bounds array is shown.
  if (cons.upper_bounds != NULL)
    {
      os << "Scaled upper bounds: (";
      for (i = 0; i < (max_dimension - 1); i++)
	os << cons.upper_bounds[i] << ", ";
      os << cons.upper_bounds[max_dimension - 1] << ")\n";
    }
  return os;
}
