algorithms-for-computing-line-estimators/LinearRegressionTool/src/main/java/de/wwwu/awolf/presenter/algorithms/advanced/LeastMedianOfSquaresEstimator.java

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.AbstractPresenter;
import de.wwwu.awolf.presenter.algorithms.Algorithm;
import de.wwwu.awolf.presenter.util.IntersectionComputer;
import de.wwwu.awolf.presenter.util.Logging;

import java.util.*;
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 LeastMedianOfSquaresEstimator implements Algorithm {
    private static final Algorithm.Type type = Type.LMS;

    private List<Line> setOfLines;
    private int n;
    private double quantileError;
    private int kPlus;
    private int kMinus;
    private Queue<Interval> 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;
    private Flow.Subscriber<? super AlgorithmData> subscriber;
    private AbstractPresenter presenter;

    /**
     * Algorithmus zum berechnen des LMS-Schätzers
     * <p>
     * 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
     */
    @Override
    public Line call() {
        this.n = setOfLines.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);

        Logging.logInfo("=== 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<Interval> 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(-10000, 10000));
        heightsigmaMin = Double.MAX_VALUE;

        //(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...
                    Collection<Point> tmpIntersections = IntersectionComputer.getInstance().compute(setOfLines, interval.getLower(), interval.getUpper());
                    boolean found = false;
                    for (Point tmpIntersection : tmpIntersections) {
                        if (tmpIntersection.getX() > interval.getLower()
                                && tmpIntersection.getX() < interval.getUpper()) {
                            intersectionsPoint = tmpIntersection.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();
        Logging.logInfo("=== E N D - L M S === " + ((end - start) / 1000));

        return pepareResult();
    }

    @Override
    public void setInput(Set<Line> lines) {
        this.setOfLines = new LinkedList<>(lines);
        setN(this.setOfLines.size());
    }

    @Override
    public Type getType() {
        return type;
    }

    @Override
    public void setPresenter(AbstractPresenter presenter) {
        this.presenter = presenter;
        subscribe(presenter);
    }


    /**
     * Liefert die Anzahl der Schnittpunkte in einem Intervall
     *
     * @param interval Intervall
     * @return Anzahl der Schnittpunkte
     */
    public int countInversions(Interval interval) {
        return IntersectionComputer.getInstance().compute(setOfLines, interval.getLower(), interval.getUpper()).size();
    }


    /**
     * 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
        List<Point> xQueue = new ArrayList<>();
        Collection<Point> points = IntersectionComputer.getInstance().compute(setOfLines, interval.getLower(), interval.getUpper());
        for (Point point : points) {
            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;
        List<Double> 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.
        List<Double> umaxList;
        List<Double> uminList;

        //y koordinaten der Schnittpunkte
        List<Line> lines = new ArrayList<>();
        for (Line p : setOfLines) {
            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 List<Double> getEjValues(double u) {

        List<Double> ret = new ArrayList<>();

        for (Line p : setOfLines) {
            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
        List<Double> intersections = new LinkedList<>();
        for (Line line : setOfLines) {
            intersections.add((px.getX() * line.getM()) + line.getB());
        }
        if (intersections.size() >= kMinusValue) {
            Collections.sort(intersections);
            double height = Math.abs(intersections.get(0) - intersections.get(kMinusValue - 1));
            return new Double[]{height, intersections.get(0), intersections.get(kMinusValue - 1)};
        } else {
            return null;
        }

    }

    private Line pepareResult() {
        if (this.subscriber != null) {
            double m = ((getSigmaMin().getX2() + getSigmaMin().getX1()) * 0.5);
            double b = (getSigmaMin().getY2() + getSigmaMin().getY1()) * -0.5;

            slope = m;
            yInterception = b;

            AlgorithmData data = new AlgorithmData();
            data.setAlgorithmType(getType());
            data.setType(SubscriberType.ALGORITHM);
            data.setLineData(new Line(m, b));
            this.subscriber.onNext(data);
        } else {
            double m = (getSigmaMin().getX2() + getSigmaMin().getX1()) * 0.5;
            double b = (getSigmaMin().getY2() + getSigmaMin().getY1()) * -0.5;

            slope = m;
            yInterception = b;
        }

        return new Line(getSlope(), getYInterception());
    }

    /**
     * @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;
    }

    @Override
    public void subscribe(Flow.Subscriber<? super Data> subscriber) {
        this.subscriber = subscriber;
    }

}