/******************************************************************************
* 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/                                      *
******************************************************************************/

// filter.cpp -- Filter class methods. The filter is used when there are
// general constraints in the optimization problem. It stores AugmentedPoint
// objects according to a rule of dominance.

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

// Private methods.
// ****************************************************************************

// This method does two things: First it finds the place where 'ap' could be
// inserted in the filter list; then if 'ap' meets certain conditions, it is
// inserted by calling the 'store(...)' method.
// - 'ap' is the Augmented Point to (possibly) insert in the filter list.
// - 'Fx' is the truth function value of 'ap'.
// - 'Hx' is the constraint violation function value of 'ap'.
int Filter::insert(const AugmentedPoint & ap, const double Fx, const double Hx)
{
  // The variable 'pos' will indicate the insertion position of 'ap' in the filter list.
  int pos = 0;
  // The variable 'success' can have these values:
  // - 0: Failure: 'ap' wasn't inserted in the filter list.
  // - 1: Partial success: 'ap' was inserted somewhere to the north-west of the
  //      infeasible incumbent.
  // - 2: Success: 'ap' was inserted somewhere to the south-west of the
  //      infeasible incumbent.
	
  // If the list is empty, the point is inserted.
  if (filter_list.empty())
		{
			if (display_factor >= 5)
				cout << "Insertion at the start: the filter list was empty.\n";
			store(ap);
			success = 2;
		}
  else
		{
			// 'current_Fx' is the Fx value for the current point in the filter list.
			// 'current_Hx' is the Hx value for the current point in the filter list.
			// 'previous_Fx' is the Fx value for the previous point in the filter list.
			double current_Fx, current_Hx, previous_Fx = INFINIT;
			// The stopping criterion is handled inside the loop.
			for (current = filter_list.begin(); current <= filter_list.end(); current++)
				{
					// Particular case: we're at the end of the filter list.
					if (current == filter_list.end())
						{
							// We check only the previous point's truth value.
							if (Fx < previous_Fx)
								{
									if (display_factor >= 5)
										cout << "Insertion at the end of the filter list.\n";
									store(ap);
									break;
								}
						}
					else   // General case.
						{
							// The current point's function values are captured.
							current_Fx = (*current)->getFx();
							current_Hx = (*current)->getHx();
							// 'Hx' is compared to the current point's constraint violation function value:
							// 1- If 'Hx' is greater than 'current_Hx', nothing is done, we go back to the
							//    top of the loop.
							// 2- If 'Hx' is inferior to 'current_Hx', 'Fx' is compared to the previous
							//    point's truth function value. If 'Fx' is inferior it's not dominated and
							//    can be inserted.
							// 3- If 'Hx' is equal to 'current_Hx', 'Fx' is compared to the current point's
							//    truth function value. If 'Fx' is inferior it's not dominated and can be inserted.
							if (((Hx < current_Hx) && (Fx < previous_Fx)) ||
									((Hx == current_Hx) && (Fx < current_Fx)))
								{
									// 'ap' is inserted in the filter list.
									if (display_factor >= 5)
										cout << "Insertion in the filter list, at position " << pos << ".\n";
									store(ap, pos);
									// We must check if the points in the rest of the list are dominated
									// (they have a superior truth function value) by the newly inserted point.
									checkDominance(Fx, pos);
									break;
								}
							// We update 'previous_Fx'.
							previous_Fx = current_Fx;
						}
					// The counter 'pos' is incremented.
					pos++;
				}
			
			// The return value is determined here.
			if (((Hx < infeasible_Hx) && (Fx <= infeasible_Fx)) ||
					((Hx == infeasible_Hx) && (Fx < infeasible_Fx)))
				success = 2;
			else
				if ((Hx < infeasible_Hx) && (Fx > infeasible_Fx))
					success = 1;
				else
					success = 0;
		}
	// If we have a new infeasible incumbent, we have to update it.
	//if (success == 2)
		//{
			//infeasible_Fx = Fx;
		//	infeasible_Hx = Hx;
		//}
  return success;
}

// This method is called after an augmented point is inserted in the filter. If
// the points that follow the inserted point have an higher truth value, they
// are dominated and thus erased from the list. */
void Filter::checkDominance(const double Fx, const int pos)
{
  int i;
  if (display_factor >= 5)
    cout << "Checking dominance, Fx = " << Fx << ", pos = " << pos << "\n";
  current = filter_list.begin();
  // The 'current' iterator is advanced to the point that follows the new inserted point.
  for (i = 0; i < pos + 1; i++)
    current++;
  // We check all the points following the inserted point, as long as we're not
  // at the end of the list (obviously) and their truth values are higher or
  // equal to 'Fx' (the inserted point's truth value, received as parameter).
  while ((current != filter_list.end()) && ((*current)->getFx() >= Fx))
		{
			if (display_factor >= 5)
				cout << "Deleting " << **current;
			// The pointed object is deleted.
			delete (*current);
			// The element is erased. The 'erase' function returns an iterator to the element
			// that had followed the erased one, that's how 'current' advances to the next element.
			current = filter_list.erase(current);
		}
}

// This method deletes the list of points from memory.
void Filter::cleanList()
{
  if (!filter_list.empty())
		{
			// The AugmentedPoint objects in 'filter_list' are destroyed.
			for (current = filter_list.begin(); current != filter_list.end(); current++)
				delete *current;
			// The memory for the pointers to AugmentedPoint objects is freed.
			filter_list.clear();
		}
}

// This method stores the AugmentedPoint 'ap' into the filter list.
void Filter::store(const AugmentedPoint & ap, const int pos)
{
  // A copy of 'ap' is made with the copy constructor.
  AugmentedPoint * temp = new AugmentedPoint(ap);
	
  // The augmented point is stored at the end of the vector.
  if (pos == -1)
		filter_list.push_back(temp);
  else
		{
			(*current) = filter_list.at(pos);
			// The copy is inserted into 'filter_list', before the 'current' iterator.
			filter_list.insert(current, temp);
		}
}

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

// Default constructor.
Filter::Filter(const Description * const des, const Parameters * const par)
: description(des), parameters(par)
{
  fx_decrease = 0.0, infeasible_Fx = INFINIT, infeasible_Hx = INFINIT;
  display_factor = success = 0;
	
  // There's no feasible incumbent at the start.
  feasible_incumbent = NULL;
  current = filter_list.begin();
}

// Destructor.
Filter::~Filter()
{
  delete feasible_incumbent;
  cleanList();
}

// This method returns the feasible incumbent (AugmentedPoint object).
const AugmentedPoint * Filter::getFeasibleIncumbent()
{
  return feasible_incumbent;
}

// This method returns the infeasible incumbent (AugmentedPoint object).
const AugmentedPoint * Filter::getInfeasibleIncumbent()
{
  if (display_factor >= 7)
	{
		cout << "\nIn 'Filter::getInfeasibleIncumbent()'\n";
		cout << "*************************************\n";
		cout << "II's f(x) = " << infeasible_Fx << ", II's h(x) = ";
		cout << infeasible_Hx << endl << endl;
	}
	
  // If the filter list is empty a null value is returned.
  if (filter_list.empty())
		return 0;
  else   // The infeasible incumbent will be returned.
		{
			int i = 0;
			std::vector<AugmentedPoint *>::iterator next = filter_list.begin();
		
			// The infeasible incumbent is found in the filter list.
			for (current = filter_list.begin(); current != filter_list.end(); current++)
				{
					next++;
			
					if (display_factor >= 7)
						{
							if (next != filter_list.end())
								{
									cout << i << ": f(x) = " << (*next)->getFx();
									cout << ", h(x) = " << (*next)->getHx() << endl;
								}
							else
								cout << "End of list reached!\n";
						}
			
					if ((next == filter_list.end()) || ((*next)->getHx() > infeasible_Hx))
					//if ( (next == filter_list.end()) || ((*next)->getHx() > infeasible_Hx) ||
						//	 (((*next)->getHx() == infeasible_Hx ) && ((*next)->getFx() <= infeasible_Fx)) )
						// We have found the infeasible incumbent.
						break;
					i++;
				}
		
			if (display_factor >= 7)
				{
					cout << endl << "Current's f(x) = " << (*current)->getFx();
					cout << ", current's h(x) = " << (*current)->getHx() << endl;
				}
		
			return (*current);
		}
}

// This method returns true if a feasible incumbent exists, false otherwise.
bool Filter::feasible() const
{
  // If 'feasible_incumbent' exists, it means that the starting point was feasible.
	if (feasible_incumbent != NULL)
		return true;
	else
		return false;
}

// This method finds out if the primary poll center is the feasible, or the infeasible incumbent.
bool Filter::isPrimaryFeasible() const
{	
  // If there are no general constraints in the problem...
  if ((description->getGeneralConstraintsNumber() == 0) || (filter_list.empty()) ||
			// ... or the objective function value of the infeasible incumbent is greater or equal to ...
			((feasible_incumbent != NULL) && (infeasible_Fx >=
			// ... (the objective function value of the feasible incumbent minus the poll center trigger).
			 (feasible_incumbent->getFx() - parameters->getFrameCenterTrigger()))))
		{
			if (display_factor >= 5)
				cout << "Primary poll center is feasible incumbent.\n";
			// The feasible incumbent will be the primary poll center.
			return true;
		}
	else
		{
			if (display_factor >= 5)
				cout << "Primary poll center is infeasible incumbent.\n";
			// The infeasible_incumbent will be the primary poll center.
			return false;
		}
}

// This method returns the feasible incumbent's f(x) value.
double Filter::getTruth() const
{
  // If 'feasible_incumbent' is null, it means that no feasible point was
  // found by NOMAD up to this point. A very large number is returned.
  if (feasible_incumbent == NULL)
    return INFINIT;
  else   // f(x) of the feasible incumbent is returned.
    return feasible_incumbent->getFx();
}

// An AugmentedPoint object is received, if it meets certain critera it will be
// added to one of the filter lists.
int Filter::populate(AugmentedPoint & trial)
{
  // We get the truth and constraint violation function values of the trial point.
  double feasible_Fx = INFINIT, trial_Fx = trial.getFx(), trial_Hx = trial.getHx();
  // The variable 'success' will be returned, with these possible values:
  // - 0: The point wasn't added to the filter.
  // - 1: The point was added to the infeasible set, but it's not better than
  //      the current infeasible incumbent.
  // - 2: The point was better than the current infeasible incumbent.
  // - 3: The point was better than the current feasible incumbent.
  success = 0;
  fx_decrease = 0.0;
	
  // Test to see if the trial point is feasible.
  if (trial_Hx == 0.0)
		{
			if (display_factor >= 7)
				cout << "Trial point is feasible.\n";
			// We get the truth function value of the feasible incumbent.
			if (feasible_incumbent != NULL)
				feasible_Fx = feasible_incumbent->getFx();
		
			// Test to see if we have a new best iterate.
			if (trial_Fx < feasible_Fx)
				{
					// Success! We have a new feasible incumbent.
					success = 3;
					// If 'feasible_incumbent' is NULL, it means that no feasible point was found up to now.
					if (feasible_incumbent == NULL)
						// The feasible incumbent is created by copying the trial point.
						feasible_incumbent = new AugmentedPoint(trial);
					else
						{
							// The decrease in f(x) is calculated.
							fx_decrease = abs(feasible_Fx - trial_Fx);
							// We copy the trial point in 'incumbent'.
							*feasible_incumbent = trial;
						}
					if (display_factor >= 7)
						cout << "Updated feasible incumbent = " << *feasible_incumbent << endl;
				}
		}
  else   // The trial point is infeasible.
		{
			// The point is not inserted in the filter list in these cases:
			// 1- If 'trial_Fx' equals 'INFINIT', the point is outside the bound constraints.
			// 2- If 'Hx' equals 'INFINIT', the point's constraint violation function value is
			//    greater than 'Hmax'.
			if ((trial_Fx != INFINIT) && (trial_Hx != INFINIT))
				success = insert(trial, trial_Fx, trial_Hx);
			if ((display_factor >= 5) && (success > 0))
				showList();
		}
  return success;
}

//
void Filter::createFeasibleIncumbent(const AugmentedPoint & trial)
{
  feasible_incumbent = new AugmentedPoint(trial);
}

//
void Filter::showFeasibleIncumbent() const
{
  // No feasible point has been found yet.
  if (feasible_incumbent == NULL)
    cout << "NA (no feasible point was found)\n";
  else   // We have a feasible point.
				 // The feasible point is written to standard output.
    cout << *feasible_incumbent;
}

//
void Filter::showFeasibleIncumbent(ofstream & fout) const
{
  // No feasible point has been found yet.
  if (feasible_incumbent == NULL)
    fout << "NA (no feasible point was found)\n";
  else   // We have a feasible point.
				 // The feasible point is written to file.
    feasible_incumbent->showResults(fout, description->getScalesAdd(),
																		description->getScalesMult());
}

// Outputs Filter's list of points.
void Filter::showList()
{
  cout << "Filter list: ";
  cout << filter_list.size() << " point(s), capacity = ";
  cout << filter_list.capacity() << "\n";
  if (filter_list.empty())
    cout << "List of points empty!\n";
  else
		{
		// We go through the list of points.
		for (current = filter_list.begin(); current != filter_list.end();
				 current++)
			// Each point is displayed.
			cout << **current;
		}
}

//
void Filter::showResults(ofstream & fout, double * scales_add,
												 double * scales_mult)
{
  feasible_incumbent->showResults(fout, scales_add, scales_mult);
}

// This method updates the infeasible incumbent after a successful iteration.
void Filter::updateInfeasibleIncumbent()
{
  // If there's no infeasible incumbent...
  if (infeasible_Hx == INFINIT)
		{
			// ... and the list is not empty.
			if (filter_list.size() > 0)
				{
					if (display_factor >= 7)
						cout << "Infeasible incumbent is last point in filter list.\n";
					// The new infeasible incumbent is the last point in the list.
					current = filter_list.end() - 1;
					// These variables are updated.
					infeasible_Fx = (*current)->getFx();
					infeasible_Hx = (*current)->getHx();
				}
		}
  else   // There is an infeasible incumbent.
		{
			if (display_factor >= 5)
				cout << "Feasible incumbent exists.\n";
			// This iterator points the element before 'current'.
			std::vector<AugmentedPoint *>::iterator previous = filter_list.begin();
			current = filter_list.begin();
		
			if ((*previous)->getHx() != infeasible_Hx)
				{
					if (display_factor >= 5)
						cout << "'previous' different than II's h(x).\n";
					// 'current' goes through the list, starting from the second element.
					while (current != filter_list.end())
						{
							if (previous == filter_list.end() - 1)
								{
									if (display_factor >= 5)
										cout << "Case 1: II updated to last point in filter list.\n";
									// These variables are updated.
									infeasible_Fx = (*previous)->getFx();
									infeasible_Hx = (*previous)->getHx();
									break;
								}
							else
								{
									// The iterator 'current' is advanced.
									current++;
									if (((*current)->getHx() > infeasible_Hx) ||
											(((*current)->getHx() == infeasible_Hx) &&
											 ((*current)->getFx() == infeasible_Fx)))
										{
											if (display_factor >= 5)
												cout << "Case 2.\n";
											// These variables are updated.
											infeasible_Hx = (*previous)->getHx();
											infeasible_Fx = (*previous)->getFx();
											break;
										}
									else
										{
											if (((*current)->getHx() == infeasible_Hx) &&
													((*current)->getFx() < infeasible_Fx))
												{
													if (display_factor >= 5)
														cout << "Case 3.\n";
													// These variables are updated.
													infeasible_Hx = (*current)->getHx();
													infeasible_Fx = (*current)->getFx();
													break;
												}
											else
												{
													previous++;
												}
										}
								}
						}
				}
			else
				{
					if ((*previous)->getFx() < infeasible_Fx)
						{
							if (display_factor >= 5)
								cout << "Case 4.\n";
							// These variables are updated.
							infeasible_Hx = (*previous)->getHx();
							infeasible_Fx = (*previous)->getFx();
						}
					else
						{
						// Nothing is done.
						if (display_factor >= 5)
							cout << "'previous' equals II's h(x).\n";
						}
				}
		}
}

// This method writes the points of the filter list in a file, to be read later
// when the filter graph is updated.
void Filter::writeInGraphFile()
{
  // The filter graph file is opened.
  ofstream fout("LIBRARY/filter_graph.txt");
  // We go through the list of points.
  for (current = filter_list.begin(); current != filter_list.end(); current++)
    // Each point's H(x) and f(x) is written in the filter graph file.
    fout << (* current)->getHx() << "  " << (* current)->getFx() << "\n";
  // The filter graph file is closed.
  fout.close();
}
