#!/usr/bin/env rdmd-dev

/** Longest Common Subsequence, typically used as a base for writing diff/compare algorithms.

    Copyright: Per Nordlöw 2014-.
    License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0).
    Authors: $(WEB Per Nordlöw)

    See_Also: https://en.wikipedia.org/wiki/Longest_common_subsequence_problem
    See_Also: https://en.wikipedia.org/wiki/Diff
*/

module nxt.lcs;

/** Longest Common Subsequence (LCS).
    See_Also: http://rosettacode.org/wiki/Longest_common_subsequence#Recursive_version
*/
T[] lcsR(T)(in T[] a,
            in T[] b) pure nothrow
{
    import std.range.primitives : empty;
    if (a.empty || b.empty)
    {
        return null;            // undefined
    }
    if (a[0] == b[0])
    {
        return a[0] ~ lcsR(a[1 .. $],
                           b[1 .. $]);
    }
    const longest = (T[] x, T[] y) => x.length > y.length ? x : y;
    return longest(lcsR(a, b[1 .. $]),
                   lcsR(a[1 .. $], b));
}

/** Longest Common Subsequence (LCS).
    Faster Dynamic Programming Version.
    See_Also: http://rosettacode.org/wiki/Longest_common_subsequence#Faster_dynamic_programming_version
*/
T[] lcsDP(T)(in T[] a,
             in T[] b) pure /* nothrow */
{
    import std.algorithm.comparison : max;

    auto L = new uint[][](a.length + 1, b.length + 1);

    foreach (immutable i; 0 .. a.length)
    {
        foreach (immutable j; 0 .. b.length)
        {
            L[i + 1][j + 1] = (a[i] == b[j]) ? (1 + L[i][j]) :
                max(L[i + 1][j], L[i][j + 1]);
        }
    }

    import core.internal.traits : Unqual;
    Unqual!T[] result;

    for (auto i = a.length, j = b.length; i > 0 && j > 0; )
    {
        if (a[i - 1] == b[j - 1])
        {
            result ~= a[i - 1];
            i--;
            j--;
        }
        else
        {
            if (L[i][j - 1] < L[i - 1][j])
                i--;
            else
                j--;
        }
    }

    import std.algorithm : reverse;
    result.reverse();           // not nothrow
    return result;
}

/** Get LCS Lengths. */
uint[] lcsLengths(R)(R xs,
                     R ys) pure nothrow
{
    import std.algorithm.comparison : max;

    auto prev = new typeof(return)(1 + ys.length);
    auto curr = new typeof(return)(1 + ys.length);

    foreach (immutable x; xs)
    {
        import std.algorithm: swap;
        swap(curr, prev);
        size_t i = 0;
        foreach (immutable y; ys)
        {
            curr[i + 1] = (x == y) ? prev[i] + 1 : max(curr[i], prev[i + 1]);
            i++;
        }
    }

    return curr;
}

void lcsDo(T)(in T[] xs,
              in T[] ys,
              bool[] xsInLCS,
              in size_t idx = 0) pure nothrow
{
    immutable nx = xs.length;
    immutable ny = ys.length;

    if (nx == 0)
        return;

    if (nx == 1)
    {
        import std.algorithm: canFind;
        if (ys.canFind(xs[0]))
            xsInLCS[idx] = true;
    }
    else
    {
        immutable mid = nx / 2;
        const xb = xs[0.. mid];
        const xe = xs[mid .. $];
        immutable ll_b = lcsLengths(xb, ys);

        import std.range: retro;
        const ll_e = lcsLengths(xe.retro,
                                ys.retro); // retro is slow with DMD

        // import std.algorithm.comparison : max;
        //immutable k = iota(ny + 1)
        //              .reduce!(max!(j => ll_b[j] + ll_e[ny - j]));
        import std.range: iota;
        import std.algorithm: minPos;
        import std.typecons: tuple;
        immutable k = iota(ny + 1).minPos!((i, j) => tuple(ll_b[i] + ll_e[ny - i]) >
                                           tuple(ll_b[j] + ll_e[ny - j]))[0];

        lcsDo(xb, ys[0 .. k], xsInLCS, idx);
        lcsDo(xe, ys[k .. $], xsInLCS, idx + mid);
    }
}

/** Longest Common Subsequence (LCS) using Hirschberg.
    Linear-Space Faster Dynamic Programming Version.

    Time Complexity: O(m*n)
    Space Complexity: O(min(m,n))

    To speed up this code on DMD remove the memory allocations from $(D lcsLengths), and
    do not use the $(D retro) range (replace it with $(D foreach_reverse))

    See_Also: https://en.wikipedia.org/wiki/Hirschberg%27s_algorithm
    See_Also: http://rosettacode.org/wiki/Longest_common_subsequence#Hirschberg_algorithm_version
*/
const(T)[] lcs(T)(in T[] xs,
                  in T[] ys) pure /*nothrow*/
{
    auto xsInLCS = new bool[xs.length];
    lcsDo(xs, ys, xsInLCS);
    import std.range : zip;
    import std.algorithm.iteration : filter, map;
    import std.array : array;
    return zip(xs, xsInLCS).filter!q{ a[1] }.map!q{ a[0] }.array; // Not nothrow.
}

string lcs(in string a,
           in string b) pure /*nothrow*/
{
    import std.string : representation, assumeUTF;
    return lcs(a.representation,
               b.representation).assumeUTF;
}

unittest
{
    immutable x = "thisisatest";
    immutable y = "testing123testing";
    immutable z = "tsitest";
    assert(z == lcsR(x, y));
    assert(z == lcsDP(x, y));
    assert(z == lcs(x, y));
    assert("" == lcs("", ""));
    assert("" == lcs(null, null));
}

unittest
{
    immutable x = [1, 2, 3];
    immutable y = [4, 5, 6];
    immutable z = [];
    assert(z == lcsR(x, y));
    assert(z == lcsDP(x, y));
    assert(z == lcs(x, y));
}

unittest
{
    immutable x = [1, 2, 3];
    immutable y = [2, 3, 4];
    immutable z = [2, 3];
    assert(z == lcsR(x, y));
    assert(z == lcsDP(x, y));
    assert(z == lcs(x, y));
}

unittest
{
    alias T = int;
    size_t n = 300;
    auto x = new T[n];
    auto y = new T[n];
    assert(lcs(x, y).length == n);
}

version(benchmark) unittest
{
    import std.algorithm : levenshteinDistance, levenshteinDistanceAndPath;
    import nxt.random_ex: randInPlaceBlockwise;

    alias T = int;
    size_t n = 300;
    auto x = new T[n];
    auto y = new T[n];

    x.randInPlaceBlockwise;
    y.randInPlaceBlockwise;

    void bLCS() { const d = lcs(x, y); }
    void bLD() { const d = levenshteinDistance(x, y); }
    void bLDAP() { const dp = levenshteinDistanceAndPath(x, y); }

    import std.datetime: benchmark;
    import std.algorithm: map;
    import std.array: array;

    auto r = benchmark!(bLCS, bLD, bLDAP)(1);
    auto q = r.array.map!(a => a.msecs);

    import std.stdio: writeln;
    writeln(q, " milliseconds");
}