//==============================================================================
//	
//	Copyright (c) 2002-
//	Authors:
//	* Vojtech Forejt <vojtech.forejt@cs.ox.ac.uk> (University of Oxford)
//	
//------------------------------------------------------------------------------
//	
//	This file is part of PRISM.
//	
//	PRISM 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 2 of the License, or
//	(at your option) any later version.
//	
//	PRISM 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 PRISM; if not, write to the Free Software Foundation,
//	Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//	
//==============================================================================

package prism;

import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

/**
 * This class contains auxiliary static methods used in multi-objective verification.
 */
public class MultiObjUtils
{
	/**
	 * Returns a point that represents a hyperplane that separates {@code point}
	 * from {@code computedPoints}. An LP solver is used to get the hyperplane, but it
	 * is expected that the number of points  in {@code computedPoints} is relatively small
	 * and so this approach should be efficient
	 * 
	 * @param point A point which should be at one direction from the separating hyperplane.
	 * @param computedPoints Set of points which should be at the other direction from the separating hyperplane.
	 * @param tolerance
	 * @return A vector orthogonal to the computed separating hyperplane.
	 * @throws PrismException When the LP solver throws an exception or returns an unexpected result, an exception with a related message is thrown.
	 */
	public static Point getWeights(Point point, List<Point> computedPoints) throws PrismException {
		int dim = point.getDimension();
		if (computedPoints == null || computedPoints.size() == 0) {
			boolean nonZero = false;
			for(int i = 0; i < dim; i++) {
				if (point.getCoord(i) != 0) {
					nonZero = true;
					break;
				}
			}			
			
			Point resPoint;
			if (nonZero) {
				resPoint = point.normalize();
			} else {
				double[] result = new double[dim];
				for(int i = 0; i < dim; i++) {
					result[i] = 1.0;
				}
				resPoint = (new Point(result)).normalize();
			}
			return resPoint;
		}
		
		int epsilonCoord = dim + 1;
		//int numVars = dim + 1;
		double[] d;
		try {
			lpsolve.LpSolve solver = lpsolve.LpSolve.makeLp(0,dim + 2);
			solver.setVerbose(lpsolve.LpSolve.CRITICAL);
			
			//for some strange reason the 0th element of an array is ignored by addConstraint
			d = new double[dim + 3];
			for (int i = 0; i< dim; i++) {
				d[i + 1] = point.getCoord(i);
			}
			
			d[dim + 2] = -1.0;
			
			//System.out.println(Arrays.toString(d));
			solver.addConstraint(d, lpsolve.LpSolve.GE, 0.0);
			
			for (Point p : computedPoints)
			{
				d = new double[dim + 3];
				for (int i = 0; i< dim; i++) {
					d[i + 1] = p.getCoord(i);
				}
				d[epsilonCoord] = 1;
				
				d[dim + 2] = -1.0;
				
				solver.addConstraint(d, lpsolve.LpSolve.LE, 0.0);
			}
			
			for(int i = 0; i < dim; i++) {
				d = new double[dim + 3];
				d[i + 1] = 1;
				
				solver.addConstraint(d, lpsolve.LpSolve.GE, 0.0);
			}
			
			d = new double[dim + 3];
			
			for(int i = 0; i < dim; i++) {
				d[i + 1] = 1;
			}

			solver.addConstraint(d, lpsolve.LpSolve.LE, 1.0);

			d = new double[dim + 3];
			d[epsilonCoord] = 1;
				
			solver.setBounds(dim + 2, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
			
			solver.setObjFn(d);
			solver.setMaxim();
			
			int r = solver.solve();
					
			if (r == lpsolve.LpSolve.INFEASIBLE) {
				System.out.println("The solution is infeasible");
				return null;
			} else if (r == lpsolve.LpSolve.UNBOUNDED) {
				//this should never happen
				throw new PrismException("The solution of a linear program is unbounded.");
			} if (r != lpsolve.LpSolve.OPTIMAL)
				//this should also never happen
				throw new PrismException("LpSolve returns an unexpected value: " + r);

			double[] resultSolver = solver.getPtrVariables();
			//solver.printLp();
			if (resultSolver[epsilonCoord-1] == 0.0) {
				//there is no _strictly_ separating hyperplane (i.e. one that doesn't touch any of the points)
				return null;
			}
			
			//solver.printLp();
		    //System.out.println("solver says:" + r);
		    //System.out.println("solver results:" + Arrays.toString(resultSolver));
				
			double[] result = new double[dim];
			//copy all coordinates and ignore auxiliary variables
			for(int i = 0; i < dim; i++)
				result[i] = resultSolver[i];
			
			Point resPoint = new Point(result).normalize();
			
			return resPoint;
		} catch (lpsolve.LpSolveException ex) {
			PrismException ex2 = new PrismException("lpsolve threw an exception: " + ex.getMessage());
			throw ex2;
		}
	}

	/**
	 * Prints an auxiliary graph file in gnuplot format. Points must be 2D.
	 * 
	 * NOTE: not very nice now, should be enhanced in the future...
	 * 
	 * @param target Current target
	 * @param points Corner points of the pareto curve
	 * @param directions Directions associated with the corner points.
	 * @param directory The directory where to store the file.
	 * @param index Counter that should be incremented with every call of this method, determines a suffix of the name of the file.
	 * @throws PrismException
	 */
	public static void printGraphFileDebug(Point target, List<Point> points, List<Point> directions, String directory, int index) throws PrismException {
		//System.out.println("Plotting " + index + ": " + points + ", " + directions);
		
		double padding = 1.2;
		double maxX = target.getCoord(0) * padding;
		double maxY = target.getCoord(1) * padding;
		for (Point p : points) {
			double x = p.getCoord(0) * padding;
			double y = p.getCoord(1) * padding;
			maxX = (maxX >= x) ? maxX : x;
			maxY = (maxY >= y) ? maxY : y;
		}
		
		try {
			FileWriter fw;
			if (target != null) {
				fw = new FileWriter(new File(directory, "dataT" + index + ".dat"));
				fw.write(target.getCoord(0) + " " + target.getCoord(1));
				fw.close();
			}
			
			fw = new FileWriter(new File(directory,"dataP" + index + ".dat"));
			for (Point p : points) {
				fw.write(p.getCoord(0) + " " + p.getCoord(1) + "\n");
			}
			fw.close();
			
			fw = new FileWriter(new File(directory, "plot" + index + ".plot"));
			fw.write("set parametric\n");
			fw.write("set trange [0:" + maxX + "]\n");
			fw.write("set xrange [0:" + maxX + "]\n");
			fw.write("set yrange [0:" + maxY + "]\n");
			fw.write("set term aqua\n");
			//fw.write("set output \"output" + index + ".ps\"\n");
			fw.write("plot ");
			if (target != null)
				fw.write("\"dataT" + index + ".dat\" with points pt 9, ");
			fw.write("\"< sort dataP" + index + ".dat\"  with linespoints pt 9 lw 3");
			for (int i = 0; i < points.size(); i++) {
				Point p = points.get(i);
				Point d = directions.get(i);
				double k = (p.getCoord(0)*d.getCoord(0)) + (p.getCoord(1)*d.getCoord(1));
				if (d.getCoord(1) > 0)
					fw.write(", t,(-" + d.getCoord(0) + "*t + " + k + ")/" + d.getCoord(1) + " with lines lc 3");
				else
					fw.write(", " + p.getCoord(0) + ",t with lines lc 3");
			}
			fw.close();
		} catch (IOException ex) {
			throw new PrismException("An IOException error occured when writing graph files (exception message: " + ex.getMessage() + ").");
		}
	}

	public static void exportPareto(TileList tl, String filename) throws PrismException {
		//System.out.println("Plotting " + index + ": " + points + ", " + directions);
		try {
			FileWriter fw;
			fw = new FileWriter(filename);
			fw.write(tl.toString());
			fw.close();
				} catch (IOException ex) {
			throw new PrismException("An IOException error occured when writing graph files (exception message: " + ex.getMessage() + ").");
		}
	}
	
	/**
	 * This method takes a list of computed points together with their associated directions, where points are possibly
	 * repeated, and returns a map in which to each point from {@code computedPoints} corresponds the associated direction
	 * from {@code directions} which has the most extreme slope. Which extreme is taken depends on the {@code horiz} parameter.
	 * <p/>
	 * Note that this method considers two points as "equal" if the method {@link #isCloseTo(Point)} returns {@code true}.
	 * So the returned map will not contain two keys which are too close to each other, which should be fine w.r.t. roundoff
	 * errors, but might theoretically cause trouble if the corner points of the Pareto curve are very close to each other.
	 * <p/>
	 * It is guaranteed that if there are two points {@code p1} and {@code p2} in {@code computedPoints} such that
	 * {@code p1.isCloseTo{p2}} is {@code true} and {@code p1} occurs in {@code computedPoints} before
	 * {@code p2}, then {@code p1} will be given precedence when picking a representative for the key in the map. This is crucial
	 * when using this method together with {@link MultiObjUtils#removeDuplicities(List)}.
	 * 
	 * @param computedPoints The computed points of the Pareto curve
	 * @param directions The directions associated with the computed points.
	 * @param horiz If true, gets the one which is closest to horizontal, if false then the one which is closest to vertical.
	 * @return A map in which to each point from {@code computedPoints} is associated the most extreme direction from {@code directions}.
	 */
	protected static java.util.Map<Point, Point> fillDirections(List<Point> computedPoints, List<Point> directions, boolean horiz) {
		HashMap<Point, Point> map = new HashMap<Point, Point>();
		for(int i = 0; i < computedPoints.size(); i++) {
			Point p = computedPoints.get(i);
			Point dir = directions.get(i);
			boolean added = false;
			for (Point key : map.keySet()) {
				if (key.isCloseTo(p)) {
					double slopeNew = (dir.getCoord(1) == 0) ? Double.POSITIVE_INFINITY :  dir.getCoord(0)/dir.getCoord(1);
					double slopeOrig = (map.get(key).getCoord(1) == 0) ? Double.POSITIVE_INFINITY : map.get(key).getCoord(0)/map.get(key).getCoord(1);
					if ((horiz && slopeNew < slopeOrig) || (!horiz && slopeNew > slopeOrig)) {
						map.put(key, dir);
						break;
					}
				}
			}
				
			if (!added) {
				map.put(p, dir);
			}
		}
		
		return map;
	}

	/**
	 * This method takes a list of points and returns a sub-list in which points that are very close to
	 * each other (where {@link #isCloseTo(Point)} is used to determine what is "very close") are removed.
	 * The order of elements is preserved.
	 * 
	 * It is guaranteed that if there are two points {@code p1} and {@code p2} in {@code list} such that
	 * {@code p1.isCloseTo{p2}} is {@code true} and {@code p1} occurs in {@code list} before
	 * {@code p2}, then {@code p1} will be given precedence when picking the point to keep. This is crucial
	 * when using this method together with {@link fillDirections}.
	 * @param list
	 * @return
	 */
	protected static List<Point> removeDuplicities(List<Point> list) {
		//TODO implement faster alg if we can be bothered.
		ArrayList<Point> ret = new ArrayList<Point>();
		if (list.size() == 0)
			return ret;
		
		ret.add(list.get(0));
		for (int i = 1; i < list.size(); i++) {
			Point p = list.get(i);
			
			//check if contains
			boolean contains = false;
			for(Point p2 : ret) {
				if (p2.isCloseTo(p)) {
					contains = true;
					break;
				}
			}
			
			if (!contains)
				ret.add(p);
		}
		return ret;
	}

	/**
	 * This method takes the set of some corner points of a pareto curve {@code computedPoints} together
	 * with associated orthogonal lines {@directions} that determine current over-approximation,
	 * and returns the points which determine the over-approximation. (These points are in fact the intersecting
	 * points of the directions which are not covered by other points.)
	 * @param computedPoints Corner points of the Pareto curve
	 * @param directions Directions (orthogonal lines) associated with the corner points
	 * @return The list of points that together with {@code computedPoints} determine the over-approximation of the pareto curve
	 */
	protected static List<Point> upperBoundPoints(List<Point> computedPoints, List<Point> directions)
	{
		java.util.Map<Point, Point> directionMapHoriz = fillDirections(computedPoints, directions, true);
		java.util.Map<Point, Point> directionMapVert = fillDirections(computedPoints, directions, false);
		
		PermutedLexicographicComparator plc = new PermutedLexicographicComparator(new int[] {0,1}, new boolean[] {true, false});
		List<Point> sortedLP = removeDuplicities(computedPoints);
		Collections.sort(sortedLP, plc);
		
		//System.out.println("LP: " + sortedLP);
		
		ArrayList<Point> ret = new ArrayList<Point>();
		for (int i =0; i < sortedLP.size() - 1; i++) {			
			Point p1 = sortedLP.get(i);
			Point p2 = sortedLP.get(i+1);
			Point d1 = directionMapVert.get(p1);
			Point d2 = directionMapHoriz.get(p2);
			
			//System.out.println("p1: " + p1 + ", p2: " + p2);
			Point p;
			if (d2.getCoord(1) > Point.SMALL_NUMBER) {		
				double k1 = p1.getCoord(0)*d1.getCoord(0) + p1.getCoord(1)* d1.getCoord(1);
				double k2 = p2.getCoord(0)*d2.getCoord(0) + p2.getCoord(1)* d2.getCoord(1);
				
				double[][] eqns = {{d1.getCoord(0), d1.getCoord(1)}, {d2.getCoord(0), d2.getCoord(1)}};
				double[] coords =MultiObjUtils.solveEqns(eqns, new double[] {k1,k2});
				
				if (coords == null) //HACK this is if two numbers are too close
					continue;
				
				//System.out.println("returned: " + Arrays.toString(coords));
				
				p = new Point(coords);
			} else if (d1.getCoord(0) < Point.SMALL_NUMBER && d2.getCoord(1) < Point.SMALL_NUMBER) {
				p = new Point(new double [] {p2.getCoord(0), p1.getCoord(1)});
				//System.out.println("returned': " + p);
			} else {
				double x = p2.getCoord(0);
				double k1 = p1.getCoord(0)*d1.getCoord(0) + p1.getCoord(1)* d1.getCoord(1);
				//System.out.println("x: " + x + ", k1: " + k1 + ", d1[0]: " + d1.getCoord(0) + ", d1[1]: " + d1.getCoord(1));
				double y = (k1 - d1.getCoord(0)*x)/d1.getCoord(1);
				p = new Point(new double[] {x,y});
				//System.out.println("returned'': " + p.toString());
			}
			
			//check if this point is in computedPoints
			boolean isCornerPoint = false;
			for (Point pc : computedPoints)
				if (pc.isCloseTo(p)) {
					isCornerPoint = true;
					break;
				}
					
			if (!isCornerPoint)
				ret.add(p);
		}
		return ret;
	}

	/**
	 * Returns new point for computation of two dimensional pareto curve. More formally,
	 * this method takes the points on the pareto curve computed so far, together with
	 * the directions (i.e. orthogonal hyperplanes) corresponding to the points, and
	 * returns a point which can still be in a pareto curve, and which is most distant
	 * from the points which are already known to be on the pareto curve. Here, most distant
	 * means that a point (x,y) is taken such that for every point (a,b) from the pareto
	 * curve the number max{x-a,y-b) is maximized.  
	 * 
	 * @param computedPoints The list of already computed corner points of the pareto curve.
	 * @param directions The list of directions associated with the computed corner points.
	 * @param errorTolerance The required precision of the pareto curve being approximated
	 * @return The point, or {@code null} if no point has a distance greater than {@code errorTolerance}
	 */
	protected static Point getNewTarget(List<Point> computedPoints, List<Point> directions, double errorTolerance, double[] maxValues)
	{
		//if we have less than 2 computed points, we take the most extreme directions
		if (computedPoints.size() == 0)
			return new Point(new double[] {maxValues[0], 0});
		if (computedPoints.size() == 1)
			return new Point(new double[] {0, maxValues[1]});
			//return new Point(new double[] {0, computedPoints.get(0).getCoord(1) + 1.0});
		
		//sortedLP is the set of computed corner points (i.e. their convex closure is an under-approx. of solution)
		//sortedUP is the set of points that toget (i.e. together with sortedLP, their convex closure is an over-approx. of solution)
		
		PermutedLexicographicComparator plc = new PermutedLexicographicComparator(new int[] {0,1}, new boolean[] {true, false});
		List<Point> sortedLP = removeDuplicities(computedPoints);
		Collections.sort(sortedLP, plc);
		System.out.println("sortedLP" + sortedLP);
		
		List<Point> sortedUP = upperBoundPoints(computedPoints, directions);
		Collections.sort(sortedUP, plc);
		System.out.println("sortedUP: " + sortedUP);
		
		//if (sortedUP.size() == 0)
		//	return null;
		int iLP = 0, iUP = 0;
		double maxVal = 0.0;
		int maxUP = -1;
		while (iLP < sortedLP.size() - 1 && iUP < sortedUP.size()) {
			if(sortedUP.get(iUP).getCoord(0) > sortedLP.get(iLP + 1).getCoord(0)) {
				iLP++;
			} else if(sortedUP.get(iUP).getCoord(0) <= sortedLP.get(iLP + 1).getCoord(0)) {
				double a1  = sortedLP.get(iLP).getCoord(0);
				double a2  = sortedLP.get(iLP).getCoord(1);
				double b1  = sortedLP.get(iLP + 1).getCoord(0);
				double b2  = sortedLP.get(iLP + 1).getCoord(1);
				
				//distance on the x axis
				double x2  = sortedUP.get(iUP).getCoord(1);
				double x1 = ((a1-b1)*(x2-b2))/(a2-b2) + b1;
				double d = Math.abs((sortedUP.get(iUP).getCoord(0)-x1)/x1);
				//double d = (sortedUP.get(iUP).getCoord(0)-x1);
				//System.out.println("a1: " + a1 + "  a2: " + a2 + "  b1: " + b1 + "  b2: " + b2 + "  x1: " + x1 + "  x2: " + x2 + "  d:" + d);
				if (maxVal < d) {
					maxVal = d;
					maxUP = iUP;
					//System.out.println("val: " + maxVal + ", point: " + sortedUP.get(iUP));
				}
				
				//distance on the y axis
				x1 = sortedUP.get(iUP).getCoord(0);
				x2 = ((a2-b2)*(x1-b1))/(a1-b1) + b2;
				d = Math.abs((sortedUP.get(iUP).getCoord(1)-x2)/x2);
				//d = (sortedUP.get(iUP).getCoord(1)-x2);
				if (maxVal < d) {
					maxVal = d;
					maxUP = iUP;
				}
				
				//distance along the y axis
				
				iUP++;
			} else {
				break;
			}
		}
		
		if (maxUP != -1 && maxVal > errorTolerance) {
			//System.out.println("distance:" + maxVal);
			return sortedUP.get(maxUP);
		} else {
			return null;
		}
	}

	/**
	 * Uses Gaussian elimination to solve equations given by {@code eqns} and {@code b}.
	 * The equations are given by eqns.x=b, where eqns is a matrix, and b is a column vector.
	 * 
	 * @return Array containing values for each variable.
	 */
	protected static double[] solveEqns(double[][] eqns, double[] b)
	{
		//System.out.println("Solving eqns");
		//using Gaussian elimination
	
		//System.out.println("y " + Arrays.deepToString(eqns) + " " + Arrays.toString(b));
		
		int numVars = eqns[0].length;
		for (int i = 0; i < numVars; i++) {
			//find pivot (the greatest number) for i-th column
			int index = i;
			for (int j = i + 1; j < numVars; j++) {
				if (eqns[j][i] > eqns[index][i])
					index = j;
			}
			
			double[] row = eqns[i];
			eqns[i] = eqns[index];
			eqns[index] = row;
			double right = b[i];
			b[i] = b[index];
			b[index] = right;
			
			
			if (Math.abs(eqns[i][i]) > Point.SMALL_NUMBER)
			{
				for (int j = i+1; j < eqns.length; j++)
				{
					double factor = eqns[j][i] / eqns[i][i];
					for (int k = i; k < numVars; k++) {
						eqns[j][k] -= factor * eqns[i][k];
					}
					b[j] -= factor * b[i];
				}
			}
			
			//System.out.println("z " + Arrays.deepToString(eqns) + " " + Arrays.toString(b));
		}		
		
		//System.out.println("x " + Arrays.deepToString(eqns) + " " + Arrays.toString(b));
		
		for(int i = numVars; i < eqns.length; i++)
		{
			if (Math.abs(eqns[i][numVars-1]) > Point.SMALL_NUMBER || Math.abs(b[i]) > Point.SMALL_NUMBER) {
				System.out.println("m " + Arrays.deepToString(eqns) + " " + Arrays.toString(b));
				return null;
			}
		}
		
		double[] res = new double[eqns[0].length];
		for (int i = eqns[0].length - 1; i >= 0; i--)
		{
			double sum = 0.0;
			for (int j = i + 1; j < eqns[0].length; j++) {
				sum += eqns[i][j]*res[j];
			}
			
			if (eqns[i][i] < Point.SMALL_NUMBER && Math.abs(b[i] - sum) < Point.SMALL_NUMBER)
			{
				res[i] = 0.0;
			}
			/*
			else if (eqns[i][i] < Point.SMALL_NUMBER && Math.abs(b[i] - sum) > Point.SMALL_NUMBER)
			{
				System.out.println("n " + Arrays.deepToString(eqns) + " " + Arrays.toString(b));
				return null;
			}*/
			else 
			{
				res[i] = (b[i] - sum ) / eqns[i][i];
			}
		}
		return res;
	}

/*	protected static double[] solveEqnsUsingLP(double[][] eqns, double[] b)
	{
		int dim = eqns[0].length;
		//lpsolve.LpSolve solver = lpsolve.LpSolve.makeLp(0,dim + 1);
		solver.setVerbose(lpsolve.LpSolve.CRITICAL);
		
		//for some strange reason the 0th element of an array is ignored by addConstraint
		
		for (int i = 0; i < eqns.length; i++) {
			double[] d = new double[dim + 1];
			for (int j = 0; j< dim; j++) {
				d[j + 1] = eqns[i][j];
			}
			//solver.addConstraint(d, lpsolve.LpSolve.EQ, b[i]);
		}
		return null;
	}*/
}
