/*
D1 version of the "Closest pair problem" from:
http://rosettacode.org/wiki/Closest_pair_problem
Translation to D1 from the nice Python version.

Compute nearest pair of 2D points using two algorithms.
The first algorithm is a "brute force" comparison of every possible pair.
The second algorithm is a "divide and conquer" one based on:
www.cs.iupui.edu/~xkzou/teaching/CS580/Divide-and-conquer-fastClosestPair.ppt
*/

import std.math: abs;
import std.string: format, writefln;

// from dlibs: www.fantascienza.net/leonardo/so/libs_d.zip
import d.func: table, min, filter, sorted, xpairs, select, xrange, splitter;
import d.string: putr;
import d.random: randInt, seed;
import d.time: clock;
import d.templates: Record;
import d.extra: derecord;


T dist(T)(Record!(Record!(T,T), Record!(T,T)) pair) {
    T dx = pair.d0.d0 - pair.d1.d0;
    T dy = pair.d0.d1 - pair.d1.d1;
    return abs(dx) + abs(dy);
}


Record!(T, Record!(Record!(T,T), Record!(T,T))) closestPair(T)(Record!(T,T)[] points) {
    alias Record!(T, T) P;
    alias Record!(P, P) PP;
    alias Record!(T, PP) DistPP;

    if (points.length < 2)
        return DistPP.init;
    auto closest = min(xpairs(points), &dist!(T));
    return DistPP(dist(closest), closest);
}


Record!(T, Record!(Record!(T,T), Record!(T,T))) fastClosestPair(T)(Record!(T,T)[] points) {
    alias Record!(T, T) P;
    alias Record!(P, P) PP;
    alias Record!(T, PP) DistPP;

    static Record!(T, PP) closestPairAux(P[] xp, P[] yp) {
        int npoints = xp.length;
        if (npoints <= 3)
            return closestPair(xp);

        P[] pl = xp[0 .. npoints / 2];
        P[] pr = xp[npoints / 2 .. $];

        auto x_divider = pl[$ - 1].d0;
        P[] yl, yr;
        splitter(yp, (P p){ return p.d0 <= x_divider; }, yl, yr);

        T dl, dr = void;
        PP pair_l, pair_r = void;
        derecord(closestPairAux(pl, yl), dl, pair_l);
        derecord(closestPairAux(pr, yr), dr, pair_r);

        T dm = (dl < dr) ? dl : dr;
        PP pair_m = (dl < dr) ? pair_l : pair_r;

        // Points within dm of x_divider sorted by y coord
        auto close_y = filter((P p){ return abs(p.d0 - x_divider) < dm; }, yp);
        int num_close_y = close_y.length;
        if (num_close_y > 1) {
            // very slow, it should be lazy as in Python!
            int i, j;
            auto pairs = select(PP(close_y[i], close_y[j]),
                                i, xrange(close_y.length - 1),
                                j, xrange(i + 1, min(i + 8, num_close_y)));

            PP closest_y = min(pairs, &dist!(T));
            return dm <= dist(closest_y) ? DistPP(dm, pair_m) : DistPP(dist(closest_y), closest_y);
        } else
            return DistPP(dm, pair_m);
    }

    P[] xp = sorted(points, (P p){ return p.d0; });
    P[] yp = sorted(points, (P p){ return p.d1; });
    return closestPairAux(xp, yp);
}


void main() {
    alias int T;
    alias Record!(T, T) P;
    seed(1);
    P[] points = [P(5,9), P(9,3), P(2,0), P(8,4), P(7,4), P(9,10), P(1,9),
                  P(8,2), P(0,10), P(9,6)];
    putr(points);

    putr("    closestPair: ", closestPair(points));
    putr("fastClosestPair: ", fastClosestPair(points));
    putr();

    // show results for some randomly generated lists of points
    for (int i; i < 10; i++) {
        points = table(P(randInt(10), randInt(10)), 10);
        putr(points);
        putr("    closestPair: ", closestPair(points));
        putr("fastClosestPair: ", fastClosestPair(points));
        putr();
    }

    // time the different functions
    for (int i; i < 3; i++) {
        points = table(P(randInt(100_000), randInt(100_000)), 2_000);
        putr();

        auto t0 = clock();
        auto dist_pair_1 = closestPair(points);
        auto t1 = clock();
        putr(dist_pair_1);
        writefln("Time for closestPair %.3f", t1 - t0);

        t0 = clock();
        auto dist_pair_2 = fastClosestPair(points);
        t1 = clock();
        putr(dist_pair_2);
        writefln("Time for fastClosestPair %.3f", t1 - t0);

        if (dist_pair_1.d0 != dist_pair_2.d0)
            throw new Exception("p1 != p2");
    }
}