package de.wwwu.awolf.presenter.algorithms.advanced; import de.wwwu.awolf.model.Interval; import de.wwwu.awolf.model.Line; import de.wwwu.awolf.model.Point; import de.wwwu.awolf.model.communication.AlgorithmData; import de.wwwu.awolf.model.communication.Data; import de.wwwu.awolf.model.communication.SubscriberType; import de.wwwu.awolf.presenter.Presenter; import de.wwwu.awolf.presenter.algorithms.Algorithm; import de.wwwu.awolf.presenter.util.BinomialCoeffizient; import de.wwwu.awolf.presenter.util.FastElementSelector; import de.wwwu.awolf.presenter.util.IntersectionCounter; import de.wwwu.awolf.presenter.util.RandomSampler; import java.util.ArrayList; import java.util.LinkedHashSet; import java.util.List; import java.util.Set; import java.util.concurrent.Flow; /** * Implementierung verschiedener Algorithmen zur Berechnung von Ausgleichsgeraden. * * @Author: Armin Wolf * @Email: a_wolf28@uni-muenster.de * @Date: 28.05.2017. */ public class TheilSenEstimator implements Algorithm, Flow.Publisher { private final double POSITIV_INF = 99999.0; private final double NEGATIV_INF = -99999.0; private final double EPSILON = 0.00001; private List setOfLines; private List setOfIntersections; private List intervalIntersections; private List sampledIntersections; //wird benötigt um den y Achsenabschnitt zu Berechnen private List yCoordinates; private List xCoordinates; //Hilfsvariablen (siehe original Paper) private double j; private int jA; private int jB; private double r; private int n; private double N; private int k; //Intervall und die temporaeren Grenzen private Interval interval; private double aVariant; private double bVariant; private double slope; private double yInterception; private Flow.Subscriber subscriber; /** * Konstruktor * * @param setOfLines Liste der Geraden * @param setOfIntersections Liste der Schnittpunkte * @param presenter Presenter (Beobachter) */ public TheilSenEstimator(List setOfLines, List setOfIntersections, Presenter presenter) { this.setOfLines = new ArrayList<>(setOfLines); this.setOfIntersections = new ArrayList<>(setOfIntersections); this.intervalIntersections = new ArrayList<>(setOfIntersections); this.n = setOfLines.size(); this.sampledIntersections = new ArrayList<>(); this.yCoordinates = new ArrayList<>(); this.xCoordinates = new ArrayList<>(); this.N = BinomialCoeffizient.run(n, 2); //this.k = Integer.valueOf((int) (N * 0.5)) - 1; this.k = (int) (N / 2); interval = new Interval(NEGATIV_INF, POSITIV_INF); subscribe(presenter); } /** * Konstruktor * * @param setOfLines Liste der Geraden * @param setOfIntersections Liste der Schnittpunkte */ public TheilSenEstimator(List setOfLines, List setOfIntersections) { this(setOfLines, setOfIntersections, null); } /** * Randomisierter Algorithmus zur Berechnung des Theil-Sen Schätzers. * Algorithmus stammt aus dem Paper: * "Jiri Matousek, Randomized optimal algorithm for slope selection, * Information Processing Letters 39 (1991) 183-187 */ public void run() { //damit eine initiale Ordnung herscht //Collections.sort(intervalIntersections); r = n; while (true) { if (this.N <= n || (Math.abs(interval.getUpper() - interval.getLower())) < EPSILON) { break; } else { //Anzahl der Schnittpunkte im Intervall [-Inf, a) int numberOfIntersections = getOpenIntervalSize(NEGATIV_INF, interval.getLower(), setOfIntersections); //Randomized Interpolating Search j = (r / N) * (double) (k - numberOfIntersections); jA = (int) Math.max(1, Math.floor(j - (1.5 * Math.sqrt(r)))); jB = (int) Math.min(r, Math.floor(j + (1.5 * Math.sqrt(r)))); /* Suche nach einem passenderen und kleineren Intervall Schleife terminiert wenn die das k-te Elemnet zwischen aVariant und bVariant liegt und das Intrvall weniger als 11*N / sqrt(r) Elemente besitzt */ do { //zufällige Stichprobe sampledIntersections = RandomSampler.run(intervalIntersections, r); aVariant = FastElementSelector.randomizedSelect(sampledIntersections, jA); bVariant = FastElementSelector.randomizedSelect(sampledIntersections, jB); } while (!checkCondition()); interval.setLower(aVariant); interval.setUpper(bVariant); intervalIntersections = getOpenIntervalElements(interval.getLower(), interval.getUpper()); N = getOpenIntervalSize(interval.getLower(), interval.getUpper(), intervalIntersections); } } pepareResult(); } /** * Diese Funktion überprüft ob die Bedingung für das Interval erfüllt ist. Dabei muss der k-te * Schnittpunkt in diesem Interval enthalten sein. des weiteren soll die Anzahl der Schnittpunkte * im Interval kleiner oder gleich dem Term: (11*N)/sqrt(r) sein. * * @return Boolscher Wert ob die Bedingung erfüllt ist */ private Boolean checkCondition() { //Double kthElement = FastElementSelector.randomizedSelect(xCoordinates, k); //Boolean cond1 = (kthElement > aVariant) && (kthElement <= bVariant); int lowerCount = getIntervalSize(NEGATIV_INF, aVariant); int higherCount = getIntervalSize(NEGATIV_INF, bVariant); Boolean conda = k > lowerCount; Boolean condb = k <= higherCount; Boolean cond1 = conda && condb; Boolean cond2 = (higherCount - lowerCount) <= ((11 * N) / Math.sqrt(r)); return (cond1 && cond2) || (aVariant == bVariant); } /** * Berechne wieviele von den Schnittpunkten in dem Interval zwischen a und b * enthalten sind. * * @param a untere Grenze des Intervals * @param b obere Grenze des Intrvals * @return Anzahl der Schnittpunkte im Interval [a,b) */ public int getIntervalSize(double a, double b) { IntersectionCounter ic = new IntersectionCounter(); return ic.run(setOfLines, new Interval(a, b)); } /** * Berechne wieviele von den Schnittpunkten in dem Interval zwischen a und b * enthalten sind. *

* Inspiriert durch: * https://stackoverflow.com/questions/136474/best-way-to-pick-a-random-subset-from-a-collection * * @param a untere Grenze des Intervals * @param b obere Grenze des Intrvals * @return Anzahl der Schnittpunkte im Interval (a,b) */ public int getOpenIntervalSize(double a, double b, List set) { int counter = 0; for (int i = 0; i < set.size(); i++) { Point x = set.get(i); if (x.getX() > a && x.getX() < b) { counter++; } } return counter; } /** * Berechne wieviele von den Schnittpunkten in dem Interval zwischen a und b * enthalten sind. Zusätzlich werden diese Schnittpunkte in einer Liste festgehalten und diese werden * zurückgeliefert. * * @param a untere Grenze des Intervals * @param b obere Grenze des Intrvals * @return Liste der Schnittpunkte die im Interval (a,b) vertreten sind */ public List getOpenIntervalElements(double a, double b) { List list = new ArrayList<>(); for (int i = 0; i < intervalIntersections.size(); i++) { Point x = intervalIntersections.get(i); if ((x.getX() > a && x.getX() < b) || (Math.abs(interval.getUpper() - interval.getLower())) < EPSILON) { list.add(x); } } intervalIntersections.clear(); intervalIntersections = null; return list; } @Override public void pepareResult() { double m, x; double b, y; List resultSt = getOpenIntervalElements(interval.getLower(), interval.getUpper()); List resultAbscissas = new ArrayList<>(); for (Point p : resultSt) { resultAbscissas.add(p.getX()); } for (Point p : setOfIntersections) { yCoordinates.add(p.getY()); } double pseudoIndex = getOpenIntervalSize(NEGATIV_INF, interval.getLower(), setOfIntersections) * 1.0; m = FastElementSelector.randomizedSelect(resultAbscissas, k - pseudoIndex); Set unique = new LinkedHashSet<>(yCoordinates); yCoordinates.clear(); yCoordinates.addAll(unique); b = FastElementSelector.randomizedSelect(yCoordinates, yCoordinates.size() * 0.5) * -1; slope = m; yInterception = b; if (this.subscriber != null) { AlgorithmData data = new AlgorithmData(); data.setType(SubscriberType.TS); data.setLineData(new Line(m, b)); this.subscriber.onNext(data); } } /** * @return Steigung */ public Double getSlope() { return slope; } /** * @return y-Achsenabschnitt */ public Double getyInterception() { return yInterception; } @Override public void subscribe(Flow.Subscriber subscriber) { this.subscriber = subscriber; } }