gt: 		for(int i = numVars; i < eqns.length; i++)
context: //==============================================================================
  //	
  //	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));
  		
  