#!/usr/bin/env rdmd-dev

/** Integer Sorting Algorithms.
    Copyright: Per Nordlöw 2018-.
    License: $(WEB boost.org/LICENSE_1_0.txt, Boost License 1.0).
    Authors: $(WEB Per Nordlöw)
 */
module nxt.integer_sorting;

import std.range.primitives : isRandomAccessRange, ElementType;
import std.traits : isNumeric;
import std.meta : AliasSeq;

import nxt.bijections;

/** Radix sort of `input`.

    Note that this implementation of non-inplace radix sort only requires
    `input` to be a `BidirectionalRange` not a `RandomAccessRange`.

    Note that `input` can be a `BidirectionalRange` aswell as
    `RandomAccessRange`.

    `radixBitCount` is the number of bits in radix (digit)

    TODO make `radixBitCount` a template parameter either 8 or 16,
    ElementType.sizeof must be a multiple of radixBitCount

    TODO input[] = y[] not needed when input is mutable

    TODO Restrict fun.

    TODO Choose fastDigitDiscardal based on elementMin and elementMax (if they
    are given)

    See_Also: https://probablydance.com/2016/12/27/i-wrote-a-faster-sorting-algorithm/
    See_Also: https://github.com/skarupke/ska_sort/blob/master/ska_sort.hpp
    See_Also: http://forum.dlang.org/thread/vmytpazcusauxypkwdbn@forum.dlang.org#post-vmytpazcusauxypkwdbn:40forum.dlang.org
 */
auto radixSort(R,
               alias fun = "a",
               bool descending = false,
               bool requestDigitDiscardal = false,
               bool inPlace = false)(R input,
                                     /* ElementType!R elementMin = ElementType!(R).max, */
                                     /* ElementType!R elementMax = ElementType!(R).min */)

    @trusted
if (isRandomAccessRange!R &&
    (isNumeric!(ElementType!R)))
{
    import std.range : assumeSorted;
    import std.algorithm.sorting : isSorted; // TODO move this to radixSort when know how map less to descending
    import std.algorithm.comparison : min, max;
    import std.range.primitives : front;

    immutable n = input.length; // number of elements
    alias E = ElementType!R;
    enum elementBitCount = 8*E.sizeof; // total number of bits needed to code each element

    /* Lookup number of radix bits from sizeof `ElementType`.
       These give optimal performance on Intel Core i7.
    */
    static if (elementBitCount == 8 ||
               elementBitCount == 24)
    {
        enum radixBitCount = 8;
    }
    else static if (elementBitCount == 16 ||
                    elementBitCount == 32 ||
                    elementBitCount == 64)
    {
        enum radixBitCount = 16;
    }
    else
    {
        static assert(0, "TODO handle element type " ~ e.stringof);
    }

    // TODO activate this: subtract min from all values and then immutable uint elementBitCount = is_min(a_max) ? 8*sizeof(E) : binlog(a_max); and add it back.
    enum digitCount = elementBitCount / radixBitCount;         // number of `digitCount` in radix `radixBitCount`
    static assert(elementBitCount % radixBitCount == 0,
                  "Precision of ElementType must be evenly divisble by bit-precision of Radix.");

    enum doDigitDiscardal = requestDigitDiscardal && digitCount >= 2;

    enum radix = cast(typeof(radixBitCount))1 << radixBitCount;    // bin count
    enum mask = radix-1;                                     // radix bit mask

    alias UE = typeof(input.front.bijectToUnsigned); // get unsigned integer type of same precision as \tparam E.

    import nxt.fixed_dynamic_array : FixedDynamicArray;

    static if (inPlace) // most-significant digit (MSD) first in-place radix sort
    {
        static assert(!descending, "TODO Implement descending version");

        foreach (immutable digitOffsetReversed; 0 .. digitCount) // for each `digitOffset` (in base `radix`) starting with least significant (LSD-first)
        {
            immutable digitOffset = digitCount - 1 - digitOffsetReversed;
            immutable digitBitshift = digitOffset*radixBitCount; // digit bit shift

            // [lowOffsets[i], highOffsets[i]] will become slices into `input`
            size_t[radix] lowOffsets; // low offsets for each bin
            size_t[radix] highOffsets; // high offsets for each bin

            // calculate counts
            foreach (immutable j; 0 .. n) // for each element index `j` in `input`
            {
                immutable UE currentUnsignedValue = cast(UE)input[j].bijectToUnsigned(descending);
                immutable i = (currentUnsignedValue >> digitBitshift) & mask; // digit (index)
                ++highOffsets[i];   // increase histogram bin counter
            }

            // bin boundaries: accumulate bin counters array
            lowOffsets[0] = 0;             // first low is always zero
            foreach (immutable j; 1 .. radix) // for each successive bin counter
            {
                lowOffsets[j] = highOffsets[j - 1]; // previous roof becomes current floor
                highOffsets[j] += highOffsets[j - 1]; // accumulate bin counter
            }
            assert(highOffsets[radix - 1] == n); // should equal high offset of last bin
        }

        // /** \em unstable in-place (permutate) reorder/sort `input`
        //  * access `input`'s elements in \em reverse to \em reuse filled caches from previous forward iteration.
        //  * \see `in_place_indexed_reorder`
        //  */
        // for (int r = radix - 1; r >= 0; --r) // for each radix digit r in reverse order (cache-friendly)
        // {
        //     while (binStat[r])  // as long as elements left in r:th bucket
        //     {
        //         immutable uint i0 = binStat[r].pop_back(); // index to first element of permutation
        //         immutable E    e0 = input[i0]; // value of first/current element of permutation
        //         while (true)
        //         {
        //             immutable int rN = (e0.bijectToUnsigned(descending) >> digitBitshift) & mask; // next digit (index)
        //             if (r == rN) // if permutation cycle closed (back to same digit)
        //                 break;
        //             immutable ai = binStat[rN].pop_back(); // array index
        //             swap(input[ai], e0); // do swap
        //         }
        //         input[i0] = e0;         // complete cycle
        //     }
        // }

        // TODO copy reorder algorithm into local function that calls itself in the recursion step
        // TODO call this local function

        assert(input.isSorted!"a < b");
    }
    else                        // standard radix sort
    {
        // non-in-place requires temporary `y`. TODO we could allocate these as
        // a stack-allocated array for small arrays and gain extra speed.
        auto tempStorage = FixedDynamicArray!E.makeUninitializedOfLength(n);
        auto tempSlice = tempStorage[];

        static if (doDigitDiscardal)
        {
            UE ors  = 0;         // digits diff(xor)-or-sum
        }

        foreach (immutable digitOffset; 0 .. digitCount) // for each `digitOffset` (in base `radix`) starting with least significant (LSD-first)
        {
            immutable digitBitshift = digitOffset*radixBitCount;   // digit bit shift

            static if (doDigitDiscardal)
            {
                if (digitOffset != 0) // if first iteration already performed we have bit statistics
                {
                    if ((! ((ors >> digitBitshift) & mask))) // if bits in digit[d] are either all \em zero or
                    {
                        continue;               // no sorting is needed for this digit
                    }
                }
            }

            // calculate counts
            size_t[radix] highOffsets; // histogram buckets count and later upper-limits/walls for values in `input`
            UE previousUnsignedValue = cast(UE)input[0].bijectToUnsigned(descending);
            foreach (immutable j; 0 .. n) // for each element index `j` in `input`
            {
                immutable UE currentUnsignedValue = cast(UE)input[j].bijectToUnsigned(descending);
                static if (doDigitDiscardal)
                {
                    if (digitOffset == 0) // first iteration calculates statistics
                    {
                        ors |= previousUnsignedValue ^ currentUnsignedValue; // accumulate bit change statistics
                        // ors |= currentUnsignedValue; // accumulate bits statistics
                    }
                }
                immutable i = (currentUnsignedValue >> digitBitshift) & mask; // digit (index)
                ++highOffsets[i];              // increase histogram bin counter
                previousUnsignedValue = currentUnsignedValue;
            }

            static if (doDigitDiscardal)
            {
                if (digitOffset == 0) // if first iteration already performed we have bit statistics
                {
                    if ((! ((ors >> digitBitshift) & mask))) // if bits in digit[d] are either all \em zero or
                    {
                        continue;               // no sorting is needed for this digit
                    }
                }
            }

            // bin boundaries: accumulate bin counters array
            foreach (immutable j; 1 .. radix) // for each successive bin counter
            {
                highOffsets[j] += highOffsets[j - 1]; // accumulate bin counter
            }
            assert(highOffsets[radix - 1] == n); // should equal high offset of last bin

            // reorder. access `input`'s elements in \em reverse to \em reuse filled caches from previous forward iteration.
            // \em stable reorder from `input` to `tempSlice` using normal counting sort (see `counting_sort` above).
            enum unrollFactor = 1;
            assert((n % unrollFactor) == 0, "TODO Add reordering for remainder"); // is evenly divisible by unroll factor
            for (size_t j = n - 1; j < n; j -= unrollFactor) // for each element `j` in reverse order. when `j` wraps around `j` < `n` is no longer true
            {
                static foreach (k; 0 .. unrollFactor) // inlined (unrolled) loop
                {
                    immutable i = (input[j - k].bijectToUnsigned(descending) >> digitBitshift) & mask; // digit (index)
                    tempSlice[--highOffsets[i]] = input[j - k]; // reorder into tempSlice
                }
            }
            assert(highOffsets[0] == 0); // should equal low offset of first bin

            static if (digitCount & 1) // if odd number of digit passes
            {
                static if (__traits(compiles, input[] == tempSlice[]))
                {
                    input[] = tempSlice[]; // faster than std.algorithm.copy() because input never overlap tempSlice
                }
                else
                {
                    import std.algorithm.mutation : copy;
                    copy(tempSlice[], input[]); // TODO use memcpy
                }
            }
            else
            {
                import std.algorithm.mutation : swap;
                swap(input, tempSlice);
            }
        }
    }

    static if (descending)
    {
        return input.assumeSorted!"a > b";
    }
    else
    {
        return input.assumeSorted!"a < b";
    }
}

version = benchmark;
// version = show;

version(benchmark)
@safe unittest
{
    version(show) import std.stdio : write, writef, writeln;

    /** Test `radixSort` with element-type `E`. */
    void test(E)(int n) @safe
    {
        version(show) writef("%8-s, %10-s, ", E.stringof, n);

        import nxt.dynamic_array : Array = DynamicArray;

        import std.traits : isIntegral, isSigned, isUnsigned;
        import nxt.random_ex : randInPlace, randInPlaceWithElementRange;
        import std.algorithm.sorting : sort, isSorted;
        import std.algorithm.mutation : SwapStrategy;
        import std.algorithm.comparison : min, max, equal;
        import std.range : retro;
        import std.datetime.stopwatch : StopWatch, AutoStart;
        auto sw = StopWatch();
        immutable nMax = 5;

        // generate random
        auto a = Array!E.withLength(n);
        static if (isUnsigned!E)
        {
            // a[].randInPlaceWithElementRange(cast(E)0, cast(E)uint.max);
            a[].randInPlace();
        }
        else
        {
            a[].randInPlace();
        }
        version(show) write("original random: ", a[0 .. min(nMax, $)], ", ");

        // standard quick sort
        auto qa = a.dup;

        sw.reset;
        sw.start();
        qa[].sort!("a < b", SwapStrategy.stable)();
        sw.stop;
        immutable sortTimeUsecs = sw.peek.total!"usecs";
        version(show) write("quick sorted: ", qa[0 .. min(nMax, $)], ", ");
        assert(qa[].isSorted);

        // reverse radix sort
        {
            auto b = a.dup;
            b[].radixSort!(typeof(b[]), "a", true)();
            version(show) write("reverse radix sorted: ", b[0 .. min(nMax, $)], ", ");
            assert(b[].retro.equal(qa[]));
        }

        // standard radix sort
        {
            auto b = a.dup;

            sw.reset;
            sw.start();
            b[].radixSort!(typeof(b[]), "b", false)();
            sw.stop;
            immutable radixTime1 = sw.peek.total!"usecs";

            version(show) writef("%9-s, ", cast(real)sortTimeUsecs / radixTime1);
            assert(b[].equal(qa[]));
        }

        // standard radix sort fast-discardal
        {
            auto b = a.dup;

            sw.reset;
            sw.start();
            b[].radixSort!(typeof(b[]), "b", false, true)();
            sw.stop;
            immutable radixTime = sw.peek.total!"usecs";

            assert(b[].equal(qa[]));

            version(show)
            {
                writeln("standard radix sorted with fast-discardal: ",
                        b[0 .. min(nMax, $)]);
            }
            version(show) writef("%9-s, ", cast(real)sortTimeUsecs / radixTime);
        }

        // inplace-place radix sort
        // static if (is(E == uint))
        // {
        //     auto b = a.dup;

        //     sw.reset;
        //     sw.start();
        //     b[].radixSort!(typeof(b[]), "b", false, false, true)();
        //     sw.stop;
        //     immutable radixTime = sw.peek.usecs;

        //     assert(b[].equal(qa[]));

        //     version(show)
        //     {
        //         writeln("in-place radix sorted with fast-discardal: ",
        //                 b[0 .. min(nMax, $)]);
        //     }
        //     writef("%9-s, ", cast(real)sortTimeUsecs / radixTime);
        // }

        version(show) writeln("");
    }

    import std.meta : AliasSeq;
    immutable n = 1_00_000;
    version(show) writeln("EType, eCount, radixSort (speed-up), radixSort with fast discardal (speed-up), in-place radixSort (speed-up)");
    foreach (immutable ix, T; AliasSeq!(byte, short, int, long))
    {
        test!T(n); // test signed
        import std.traits : Unsigned;
        test!(Unsigned!T)(n); // test unsigned
    }
    test!float(n);
    test!double(n);
}

@safe unittest
{
    import std.meta : AliasSeq;

    immutable n = 10_000;

    foreach (ix, E; AliasSeq!(byte, ubyte,
                              short, ushort,
                              int, uint,
                              long, ulong,
                              float, double))
    {
        import nxt.dynamic_array : Array = DynamicArray;
        import std.algorithm.sorting : sort, isSorted;
        import std.algorithm.mutation : swap;
        import nxt.random_ex : randInPlace;

        auto a = Array!E.withLength(n);

        a[].randInPlace();
        auto b = a.dup;
        assert(a[].radixSort() == b[].sort());

        swap(a, b);
    }
}