package presenter.algorithms.advanced; import model.Interval; import model.Line; import model.Point; import presenter.Presenter; import presenter.algorithms.Algorithm; import presenter.util.IntersectionCounter; import java.util.*; /** * Implementierung verschiedener Algorithmen zur Berechnung von Ausgleichsgeraden. * * @Author: Armin Wolf * @Email: a_wolf28@uni-muenster.de * @Date: 28.05.2017. */ public class LeastMedianOfSquaresEstimator extends Observable implements Algorithm { private Presenter presenter; private LinkedList set = new LinkedList<>(); private ArrayList intersections = new ArrayList<>(); private IntersectionCounter invCounter = new IntersectionCounter(); private int n; private double quantileError; private int kPlus; private int kMinus; private PriorityQueue intervals; private Interval subSlabU1; private Interval subSlabU2; private Line sigmaMin; private double heightsigmaMin; private double intersectionsPoint; private double constant = 0.5; private double slope; private double yInterception; /** * Konstruktor * @param set Liste der Geraden * @param intersections Liste der Schnittpunkte * @param presenter Presenter (Beobachter) */ public LeastMedianOfSquaresEstimator(LinkedList set, ArrayList intersections, Presenter presenter) { this.set = set; this.intersections = intersections; n = set.size(); double qPlus = 0.5; quantileError = 0.1; double qMinus = qPlus * (1 - quantileError); kMinus = (int) Math.ceil(n * qMinus); kPlus = (int) Math.ceil(n * qPlus); this.presenter = presenter; } /** * Konstruktor * @param set Liste der Geraden * @param intersections Liste der Schnittpunkte */ public LeastMedianOfSquaresEstimator(LinkedList set, ArrayList intersections) { this(set, intersections, null); } /** * Algorithmus zum berechnen des LMS-Schätzers * * Paper: * Mount, David M, Nathan S Netanyahu, Kathleen Romanik, Ruth Silverman und Angela Y Wu * „A practical approximation algorithm for the LMS line estimator“. 2007 * Computational statistics & data Analysis 51.5, S. 2461–2486 */ public void run() { System.out.println("=== S T A R T - L M S ==="); long start; long end; start = System.currentTimeMillis(); //(2.) Let U <- (-inf, inf) be the initial active intervals... Comparator comparator = (o1, o2) -> { if (o1.getDistance() < o2.getDistance()) { return -1; } if (o1.getDistance() > o2.getDistance()) { return 1; } else { return 0; } }; intervals = new PriorityQueue<>(comparator); intervals.add(new Interval(-100000, 100000)); heightsigmaMin = Double.MAX_VALUE; ArrayList tmpIntersections = intersections; //(3.) Apply the following steps as long as the exists active intervals boolean active = true; Interval interval; while (!this.intervals.isEmpty()) { interval = this.intervals.peek(); if (interval.getActivity()) { //(a.) Select any active Interval and calc. the inversions int numberOfIntersections = countInversions(interval); //(b.) apply plane sweep if ((constant * n) >= numberOfIntersections) { sigmaMin = pseudoSweep(interval); } else { //(c.) otherwise.... // get random intersections point... Collections.shuffle(tmpIntersections, new Random()); boolean found = false; for (int i = 0; i < tmpIntersections.size(); i++) { if (tmpIntersections.get(i).getX() > interval.getLower() && tmpIntersections.get(i).getX() < interval.getUpper()) { intersectionsPoint = tmpIntersections.get(i).getX(); found = true; break; } } if (found) { splitActiveSlab(intersectionsPoint, interval); //(d.) this may update sigma min upperBound(intersectionsPoint); //(e.) for i={1,2}, call lower bound(Ui) lowerBound(subSlabU1); lowerBound(subSlabU2); if (subSlabU1.getActivity()) { this.intervals.add(subSlabU1); } if (subSlabU2.getActivity()) { this.intervals.add(subSlabU2); } } else { this.intervals.poll(); } } } else { this.intervals.remove(interval); } } end = System.currentTimeMillis(); System.out.println("Zeit: "+ ((end-start)/1000)); pepareResult(); } /** * Liefert die Anzahl der Schnittpunkte in einem Intervall * @param interval Intervall * @return Anzahl der Schnittpunkte */ public int countInversions(Interval interval) { return invCounter.run(set, interval); } /** * In der Literatur wird ein planesweep vorrausgesetzt um in einem Intervall das minimale kMinus Segment zu * identifizieren. An dieser Stelle wird eine naivere Lösung verwendet, mithilfe der Schleife werden alle * Schnittpunkte betrachtet und bei passenden x-Koordinaten werden wird die vertikale Gerade auf kMinus Segmente * untersucht. * * @param interval * @return */ public Line pseudoSweep(Interval interval) { //initialisiere die x-Queue mit den 2D Punkten und sortiere nach x-Lexikographischer Ordnung ArrayList xQueue = new ArrayList<>(); for (Point point : intersections) { if (point.getX() >= interval.getLower() && point.getX() < interval.getUpper()) { xQueue.add(point); } } Collections.sort(xQueue); Line bracelet = sigmaMin; double heightOfBracelet = heightsigmaMin; for (Point current : xQueue) { Double[] currentBracelet = calcKMinusBracelet(current, kMinus); if (currentBracelet == null) { continue; } else if (currentBracelet[0] < heightOfBracelet) { heightOfBracelet = currentBracelet[0]; bracelet = new Line(current.getX(), current.getX(), currentBracelet[1], currentBracelet[2]); } } interval.setActivity(false); return bracelet; } /** * Diese Methode spaltet den aktiven Interval an der x Koordinate point. Es werden zwei neue Slabs * erzeugt. * * @param point x Koordinate an der, der Split geschieht. */ public void splitActiveSlab(double point, Interval active) { subSlabU1 = new Interval(active.getLower() + 0.01, point); subSlabU2 = new Interval(point, active.getUpper()); this.intervals.remove(active); } /** * Zu einer vertikalen Gerade werden kMinus Segmente erzeugt und geprüft ob eine bessere Lösung * als SigmaMin vorliegt. * @param point x-Koordinate zur Konstruktion der vertikalen Gerade */ public void upperBound(double point) { double height = heightsigmaMin; double tmpHeight; ArrayList sortedLineSequence = getEjValues(point); int itnbr = ((n - kMinus) + 1); for (int i = 0; i < itnbr; i++) { tmpHeight = sortedLineSequence.get((i + kMinus) - 1) - sortedLineSequence.get(i); if (tmpHeight < height) { height = tmpHeight; } if (height < heightsigmaMin) { heightsigmaMin = height; if (sigmaMin != null) { sigmaMin.setEndPoints(point, sortedLineSequence.get(i) , point, sortedLineSequence.get((i + kMinus) - 1)); } else { sigmaMin = new Line(point, point, sortedLineSequence.get(i), sortedLineSequence.get((i + kMinus) - 1)); } } } } /** * Die Methode prüft ob das übergebene Intervall die untere Schranke einer potenziellen Lösung erfüllt. * @param pslab übergebene Intervall */ public void lowerBound(Interval pslab) { int[] alpha = new int[n]; int[] beta = new int[n]; int strictlyGreater = 0; //Teil I. ArrayList umaxList; ArrayList uminList; //y koordinaten der Schnittpunkte ArrayList lines = new ArrayList<>(); for (Line p : set) { lines.add( new Line(pslab.getLower(), pslab.getUpper(), ((pslab.getLower() * p.getM()) + p.getB()), ((pslab.getUpper() * p.getM()) + p.getB()))); } umaxList = getEjValues(pslab.getUpper()); uminList = getEjValues(pslab.getLower()); for (int i = 0; i < n; i++) { Line level = new Line(pslab.getLower(), pslab.getUpper(), uminList.get(i), umaxList.get(i)); for (Line line : lines) { if ((line.getY1() < level.getY1()) && (line.getY2() < level.getY2())) { alpha[i]++; } if ((line.getY1() > level.getY1()) && (line.getY2() > level.getY2())) { strictlyGreater++; } } beta[i] = n - (alpha[i] + strictlyGreater); strictlyGreater = 0; } //Teil II. int i = 0; double h; pslab.setActivity(false); for (int j = 0; j < n; j++) { while ((i < n && (Math.abs(beta[i] - alpha[j]) < kPlus))) { i++; } if (i >= n) { pslab.setActivity(false); break; } else { h = Math.min(Math.abs(uminList.get(j) - uminList.get(i)), Math.abs(umaxList.get(j) - umaxList.get(i))); double error = 0.01; if (((1 + error) * h) < heightsigmaMin) { pslab.setActivity(true); return; } } i = 0; } } /** * Berechnet die Schnittpunkte der Geraden und der vertikalen Gerade u. Im paper sind diese Werte * als e_j Werte bekannt. * * @param u vertikale Gerade * @return Liste der Schnittpunkte (da u bekannt werden nur die y Werte zurück gegeben) */ public ArrayList getEjValues(double u) { ArrayList ret = new ArrayList<>(); for (Line p : set) { ret.add((p.getM() * u) + p.getB()); } Collections.sort(ret); return ret; } /** * Die Funktion berechnet anhand einer vertikalen Gerade x = px das sogenannte kleinste kMinus * Bracelet. Mit anderen Worten es wird eine vertikale Teilgerade berechnet die mindestens kMinus * Geraden schneidet und dabei minimal ist. * * @param px Koordinate um die "vertikale Gerade" zu simulieren. * @return Das Array enthält höhe des Bracelet, e_j und e_(j + kMinus - 1) */ public Double[] calcKMinusBracelet(Point px, int kMinusValue) { //y Koordinaten für das kMinus brecalet LinkedList intersections = new LinkedList<>(); for (Line line : set) { intersections.add((px.getX() * line.getM()) + line.getB()); } if (intersections.size() >= kMinusValue) { Collections.sort(intersections); double height = Math.abs(intersections.get(0) - intersections.get(0 + kMinusValue - 1)); Double[] ret = {height, intersections.get(0), intersections.get(0 + kMinusValue - 1)}; return ret; } else { return null; } } @Override public void pepareResult() { if (presenter != null) { setChanged(); double m = (getSigmaMin().getX2() + getSigmaMin().getX1()) * 0.5; double b = (getSigmaMin().getY2() + getSigmaMin().getY1()) * -0.5; slope = m; yInterception = b; String[] result = {"lms", m + "", b + ""}; notifyObservers(result); } else { double m = (getSigmaMin().getX2() + getSigmaMin().getX1()) * 0.5; double b = (getSigmaMin().getY2() + getSigmaMin().getY1()) * -0.5; slope = m; yInterception = b; } } /** * @return Anzahl der Geraden */ public int getN() { return n; } /** * @param n Anzahl der Geraden */ public void setN(int n) { this.n = n; } /** * @param quantileError Parameter quantile Fehler */ public void setQuantileError(double quantileError) { this.quantileError = quantileError; } /** * @param kMinus kMinus index */ public void setkMinus(int kMinus) { this.kMinus = kMinus; } /** * @return duale Streifen Sigma_min */ public Line getSigmaMin() { return sigmaMin; } /** * @param sigmaMin duale Streifen Sigma_min */ public void setSigmaMin(Line sigmaMin) { this.sigmaMin = sigmaMin; } /** * @param heightsigmaMin höhe von Sigma_min */ public void setHeightsigmaMin(double heightsigmaMin) { this.heightsigmaMin = heightsigmaMin; } /** * @param constant Parameter Konstante */ public void setConstant(Double constant) { this.constant = constant; } /** * @return Steigung */ public Double getSlope() { return slope; } /** * @return y-Achsenabschnitt */ public Double getyInterception() { return yInterception; } }