/**Relatively low-level primitives on which to build higher-level math/stat
 * functionality.  Some are used internally, some are just things that may be
 * useful to users of this library.  This module is starting to take on the
 * appearance of a small utility library.
 *
 * Note:  In several functions in this module that return arrays, the last
 * parameter is an optional buffer for storing the return value.  If this
 * parameter is ommitted or the buffer is not large enough, one will be
 * allocated on the GC heap.
 *
 * Author:  David Simcha*/
 /*
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 ** The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */
module dstats.base;

import std.math, std.mathspecial, std.traits, std.typecons, std.algorithm,
    std.range, std.exception, std.conv, std.functional, std.typetuple,
    std.numeric;

import dstats.alloc, dstats.sort;

// Returns the number of dimensions in an array T.
package template nDims(T)
{
    static if(isArray!T)
    {
        enum nDims = 1 + nDims!(typeof(T.init[0]));
    }
    else
    {
        enum nDims = 0;
    }
}

version(unittest)
{
    static assert(nDims!(uint[]) == 1);
    static assert(nDims!(float[][]) == 2);
}

import std..string : strip;

immutable double[] logFactorialTable;

private enum size_t staticFacTableLen = 10_000;

shared static this() {
    // Allocating on heap instead of static data segment to avoid
    // false pointer GC issues.
    double[] sfTemp = new double[staticFacTableLen];
    sfTemp[0] = 0;
    for(uint i = 1; i < staticFacTableLen; i++) {
        sfTemp[i] = sfTemp[i - 1] + log(i);
    }
    logFactorialTable = assumeUnique(sfTemp);
}

version(unittest) {
    import std.stdio, std.random, std.file;
}

/**
This is the exception that is thrown on invalid arguments to a dstats function.
*/
class DstatsArgumentException : Exception {
    this(string msg, string file, int line) {
        super(msg, file, line);
    }
}

T dstatsEnforce(T, string file = __FILE__, int line = __LINE__)
(T value, lazy const(char)[] msg = null) {
    if(!value) {
        const(char)[] lazyMsg = msg;
        auto exceptMsg = (lazyMsg !is null) ? lazyMsg.idup : "Invalid argument.";
        throw new DstatsArgumentException(exceptMsg, file, line);
    }

    return value;
}

package void enforceConfidence
(string file = __FILE__, int line = __LINE__)(double conf) {
    dstatsEnforce!(bool, file, line)(conf >= 0 && conf <= 1,
        text("Confidence intervals must be between 0 and 1, not ", conf, "."));
}

/**
Tests whether T is an input range whose elements can be implicitly
converted to doubles.*/
template doubleInput(T) {
    enum doubleInput = isInputRange!(T) && is(ElementType!(T) : double);
}

/**Tests whether T is iterable and has elements of a type implicitly
 * convertible to double.*/
template doubleIterable(T) {
    static if(!isIterable!T) {
        enum doubleIterable = false;
    } else {
        enum doubleIterable = is(ForeachType!(T) : double);
    }
}

/**
Given a tuple possibly containing forward ranges, returns a tuple where
save() has been called on all forward ranges.
 */
Tuple!T saveAll(T...)(T args) {
    Tuple!T ret;

    foreach(ti, elem; args) {
        static if(isForwardRange!(typeof(elem))) {
            ret.field[ti] = elem.save;
        } else {
            ret.field[ti] = elem;
        }
    }

    return ret;
}

/**Bins data into nbin equal width bins, indexed from
 * 0 to nbin - 1, with 0 being the smallest bin, etc.
 * The values returned are the counts for each bin.
 *
 * Works with any forward range with elements implicitly convertible to double.
 */
Ret[] binCounts(Ret = uint, T)(T data, uint nbin, Ret[] buf = null)
if(isForwardRange!(T) && doubleInput!(T)) {
    dstatsEnforce(nbin > 0, "Cannot bin data into zero bins.");

    alias Unqual!(ElementType!(T)) E;
    E min = data.front, max = data.front;
    foreach(elem; data) {
        if(elem > max)
            max = elem;
        else if(elem < min)
            min = elem;
    }
    E range = max - min;

    Ret[] bins;
    if(buf.length < nbin) {
        bins = new Ret[nbin];
    } else {
        bins = buf[0..nbin];
        bins[] = 0;
    }

    foreach(elem; data) {
        // Using the truncation as a feature.
        uint whichBin = cast(uint) ((elem - min) * nbin / range);

        // Handle edge case by putting largest item in largest bin.
        if(whichBin == nbin) whichBin--;

        bins[whichBin]++;
    }

    return bins;
}

unittest {
    double[] data = [0.0, .01, .03, .05, .11, .31, .51, .71, .89, 1];
    auto res = binCounts(data, 10);
    assert(res == [4U, 1, 0, 1, 0, 1, 0, 1, 1, 1]);

    auto buf = new uint[10];
    foreach(ref elem; buf) {
        elem = uniform(0, 34534);
    }

    res = binCounts(data, 10, buf);
    assert(res == [4U, 1, 0, 1, 0, 1, 0, 1, 1, 1]);
}

/**Bins data into nbin equal width bins, indexed from
 * 0 to nbin - 1, with 0 being the smallest bin, etc.
 * The values returned are the bin index for each element.
 *
 * Default return type is ubyte, because in the dstats.infotheory,
 * entropy() and related functions specialize on ubytes, and become
 * substandially faster.  However, if you're using more than 255 bins,
 * you'll have to provide a different return type as a template parameter.*/
Ret[] bin(Ret = ubyte, T)(T data, uint nbin, Ret[] buf = null)
if(isForwardRange!(T) && doubleInput!(T) && isIntegral!(Ret)) {
    dstatsEnforce(nbin > 0, "Cannot bin data into zero bins.");
    dstatsEnforce(nbin <= (cast(uint) Ret.max + 1), "Cannot bin into " ~
        to!string(nbin) ~ " bins and store the results in a " ~
        Ret.stringof ~ ".");
    alias ElementType!(T) E;
    Unqual!(E) min = data.front, max = data.front;

    auto dminmax = data;
    dminmax.popFront;
    foreach(elem; dminmax) {
        if(elem > max)
            max = elem;
        else if(elem < min)
            min = elem;
    }
    E range = max - min;
    Ret[] bins;
    if(buf.length < data.length) {
        bins = uninitializedArray!(Ret[])(data.length);
    } else {
        bins = buf[0..data.length];
    }

    foreach(i, elem; data) {
        // Using the truncation as a feature.
        uint whichBin = cast(uint) ((elem - min) * nbin / range);

        // Handle edge case by putting largest item in largest bin.
        if(whichBin == nbin)
            whichBin--;

        bins[i] = cast(Ret) whichBin;
    }

    return bins;
}

unittest {
    auto alloc = newRegionAllocator();
    double[] data = [0.0, .01, .03, .05, .11, .31, .51, .71, .89, 1];
    auto res = bin(data, 10);
    assert(res == to!(ubyte[])([0, 0, 0, 0, 1, 3, 5, 7, 8, 9]));

    auto buf = new ubyte[20];
    foreach(ref elem; buf) {
        elem = cast(ubyte) uniform(0U, 255);
    }

    res = bin(data, 10, buf);
    assert(res == to!(ubyte[])([0, 0, 0, 0, 1, 3, 5, 7, 8, 9]));

    // Make sure this throws:
    try {
        auto foo = bin( seq(0U, 1_000U), 512);
        assert(0);
    } catch(Exception e) {
        // It's supposed to throw.
    }
}

/**Bins data into nbin equal frequency bins, indexed from
 * 0 to nbin - 1, with 0 being the smallest bin, etc.
 * The values returned are the bin index for each element.
 *
 * Default return type is ubyte, because in the dstats.infotheory,
 * entropy() and related functions specialize on ubytes, and become
 * substandially faster.  However, if you're using more than 256 bins,
 * you'll have to provide a different return type as a template parameter.*/
Ret[] frqBin(Ret = ubyte, T)(T data, uint nbin, Ret[] buf = null)
if(doubleInput!(T) && isForwardRange!(T) && hasLength!(T) && isIntegral!(Ret)) {
    dstatsEnforce(nbin > 0, "Cannot bin data into zero bins.");
    dstatsEnforce(nbin <= data.length,
        "Cannot equal frequency bin data into more than data.length bins.");
    dstatsEnforce(nbin <= (cast(ulong) Ret.max) + 1UL, "Cannot bin into " ~
        to!string(nbin) ~ " bins and store the results in a " ~
        Ret.stringof ~ ".");

    Ret[] result;
    if(buf.length < data.length) {
        result = uninitializedArray!(Ret[])(data.length);
    } else {
        result = buf[0..data.length];
    }

    auto alloc = newRegionAllocator();
    auto perm = alloc.uninitializedArray!(size_t[])(data.length);

    foreach(i, ref e; perm) {
        e = i;
    }

    static if(isRandomAccessRange!(T)) {
        bool compare(size_t lhs, size_t rhs) {
            return data[lhs] < data[rhs];
        }

        qsort!compare(perm);
    } else {
        auto dd = alloc.array(data);
        qsort(dd, perm);
        alloc.freeLast();
    }

    auto rem = data.length % nbin;
    Ret bin = 0;
    auto i = 0, frq = data.length / nbin;
    while(i < data.length) {
        foreach(j; 0..(bin < rem) ? frq + 1 : frq) {
            result[perm[i++]] = bin;
        }
        bin++;
    }
    return result;
}

unittest {
    double[] data = [5U, 1, 3, 8, 30, 10, 7];
    auto res = frqBin(data, 3);
    assert(res == to!(ubyte[])([0, 0, 0, 1, 2, 2, 1]));
    data = [3, 1, 4, 1, 5, 9, 2, 6, 5];

    auto buf = new ubyte[32];
    foreach(i, ref elem; buf) {
        elem = cast(ubyte) i;
    }

    res = frqBin(data, 4, buf);
    assert(res == to!(ubyte[])([1, 0, 1, 0, 2, 3, 0, 3, 2]));
    data = [3U, 1, 4, 1, 5, 9, 2, 6, 5, 3, 4, 8, 9, 7, 9, 2];
    res = frqBin(data, 4);
    assert(res == to!(ubyte[])([1, 0, 1, 0, 2, 3, 0, 2, 2, 1, 1, 3, 3, 2, 3, 0]));

    // Make sure this throws:
    try {
        auto foo = frqBin( seq(0U, 1_000U), 512);
        assert(0);
    } catch(Exception e) {
        // It's supposed to throw.
    }
}

/**Generates a sequence from [start..end] by increment.  Includes start,
 * excludes end.  Does so eagerly as an array.
 *
 * Examples:
 * ---
 * auto s = seq(0, 5);
 * assert(s == [0, 1, 2, 3, 4]);
 * ---
 */
CommonType!(T, U)[] seq(T, U, V = uint)(T start, U end, V increment = 1U) {
    dstatsEnforce(increment > 0, "Cannot have a seq increment <= 0.");
    dstatsEnforce(end >= start, "End must be >= start in seq.");

    alias CommonType!(T, U) R;
    auto output = uninitializedArray!(R[])
        (cast(size_t) ((end - start) / increment));

    size_t count = 0;
    for(T i = start; i < end; i += increment) {
        output[count++] = i;
    }
    return output;
}

unittest {
    auto s = seq(0, 5);
    assert(s == [0, 1, 2, 3, 4]);
}

/**Given an input array, outputs an array containing the rank from
 * [1, input.length] corresponding to each element.  Ties are dealt with by
 * averaging.  This function does not reorder the input range.
 * Return type is float[] by default, but if you are sure you have no ties,
 * ints can be used for efficiency (in which case ties will not be averaged),
 * and if you need more precision when averaging ties, you can use double or
 * real.
 *
 * Works with any input range.
 *
 * Examples:
 * ---
 * uint[] test = [3, 5, 3, 1, 2];
 * assert(rank!("a < b", float)(test) == [3.5f, 5f, 3.5f, 1f, 2f]);
 * assert(test == [3U, 5, 3, 1, 2]);
 * ---*/
Ret[] rank(alias compFun = "a < b", Ret = double, T)(T input, Ret[] buf = null)
if(isInputRange!(T) && is(typeof(input.front < input.front))) {
    auto alloc = newRegionAllocator();

    static if(!isRandomAccessRange!(T) || !hasLength!(T)) {
        return rankSort!(compFun, Ret)( alloc.array(input), buf);
    } else {

        /* It's faster to duplicate the input and then use rankSort on the
         * duplicate, but more space efficient to create an index array.  Check
         * TempAlloc to see whether we can fit everything we need on the current
         * frame.  If yes, use the faster algorithm since we're effectively
         * not saving any space by using the space-efficient algorithm.
         * If no, use the more space-efficient algorithm.
         */
        auto bytesNeeded = input.length * (ElementType!(T).sizeof + size_t.sizeof);
        if(bytesNeeded < alloc.segmentSlack) {
            return rankSort!(compFun, Ret)( alloc.array(input), buf);
        } else {
            return rankUsingIndex!(compFun, Ret)(input, buf);
        }
    }

    assert(0);
}

/* Space-efficient but slow way of computing ranks.  The extra indirection
 * caused by creating the index array wreaks havock with CPU caches, especially
 * for large arrays that don't fit entirely in cache.
 */
private Ret[] rankUsingIndex(alias compFun, Ret, T)(T input, Ret[] buf) {
    auto alloc = newRegionAllocator();
    size_t[] indices = alloc.uninitializedArray!(size_t[])(input.length);
    foreach(i, ref elem; indices) {
        elem = i;
    }

    bool compare(size_t lhs, size_t rhs) {
        alias binaryFun!(compFun) innerComp;
        return innerComp(input[lhs], input[rhs]);
    }

    dstats.sort.qsort!compare(indices);

    Ret[] ret;
    if(buf.length < indices.length) {
        ret = uninitializedArray!(Ret[])(indices.length);
    } else {
        ret = buf[0..indices.length];
    }

    foreach(i, index; indices) {
        ret[index] = i + 1;
    }

    auto myIndexed = Indexed!(T)(input, indices);

    static if(!isIntegral!Ret) {
        averageTies(myIndexed, ret, indices);
    }

    return ret;
}

private struct Indexed(T) {
    T someRange;
    size_t[] indices;

    ElementType!T opIndex(size_t index) {
        return someRange[indices[index]];
    }

    @property size_t length() {
        return indices.length;
    }
}

/**Same as rank(), but also sorts the input range.
 * The array returned will still be identical to that returned by rank(), i.e.
 * the rank of each element will correspond to the ranks of the elements in the
 * input array before sorting.
 *
 * Examples:
 * ---
 * uint[] test = [3, 5, 3, 1, 2];
 * assert(rankSort(test) == [3.5, 5, 3.5, 1.0, 2.0]);
 * assert(test == [1U, 2, 3, 4, 5]);
 * ---
 */
Ret[] rankSort(alias compFun = "a < b", Ret = double, T)(T input, Ret[] buf = null)
if(isRandomAccessRange!(T) && hasLength!(T) && is(typeof(input.front < input.front))) {
    auto alloc = newRegionAllocator();
    Ret[] ranks;
    if(buf.length < input.length) {
        ranks = uninitializedArray!(Ret[])(input.length);
    } else {
        ranks = buf[0..input.length];
    }

    size_t[] perms = alloc.uninitializedArray!(size_t[])(input.length);
    foreach(i, ref p; perms) {
        p = i;
    }

    dstats.sort.qsort!compFun(input, perms);
    foreach(i; 0..perms.length)  {
        ranks[perms[i]] = i + 1;
    }

    static if(!isIntegral!Ret) {
        averageTies(input, ranks, perms);
    }

    return ranks;
}

unittest {
    uint[] test = [3, 5, 3, 1, 2];

    float[] dummy;
    assert(rankUsingIndex!("a < b", float)(test, dummy) == [3.5f, 5f, 3.5f, 1f, 2f]);
    assert(test == [3U, 5, 3, 1, 2]);

    double[] dummy2;
    assert(rankUsingIndex!("a < b", double)(test, dummy2) == [3.5, 5, 3.5, 1, 2]);
    assert(rankSort(test) == [3.5, 5.0, 3.5, 1.0, 2.0]);
    assert(test == [1U,2,3,3,5]);

    uint[] test2 = [3,3,1,2];
    assert(rank(test2) == [3.5,3.5,1,2]);
}

private void averageTies(T, U)(T sortedInput, U[] ranks, size_t[] perms)
in {
    assert(sortedInput.length == ranks.length);
    assert(ranks.length == perms.length);
} body {
    size_t tieCount = 1;
    foreach(i; 1..ranks.length) {
        if(sortedInput[i] == sortedInput[i - 1]) {
            tieCount++;
        } else if(tieCount > 1) {
            double avg = 0;
            immutable increment = 1.0 / tieCount;

            foreach(perm; perms[i - tieCount..i]) {
                avg += ranks[perm] * increment;
            }

            foreach(perm; perms[i - tieCount..i]) {
                ranks[perm] = avg;
            }
            tieCount = 1;
        }
    }

    if(tieCount > 1) { // Handle the end.
        double avg = 0;
        immutable increment = 1.0 / tieCount;

        foreach(perm; perms[perms.length - tieCount..$]) {
            avg += ranks[perm] * increment;
        }

        foreach(perm; perms[perms.length - tieCount..$]) {
            ranks[perm] = avg;
        }
    }
}

/**Returns an associative array of counts of every element in input.
 * Works w/ any iterable.
 *
 * Examples:
 * ---
 * int[] foo = [1,2,3,1,2,4];
 * uint[int] frq = frequency(foo);
 * assert(frq.length == 4);
 * assert(frq[1] == 2);
 * assert(frq[4] == 1);
 * ---*/
uint[ForeachType!(T)] frequency(T)(T input)
if(isIterable!(T)) {
    typeof(return) output;
    foreach(i; input) {
        output[i]++;
    }
    return output;
}

unittest {
    int[] foo = [1,2,3,1,2,4];
    uint[int] frq = frequency(foo);
    assert(frq.length == 4);
    assert(frq[1] == 2);
    assert(frq[4] == 1);
}

/**Given a range of values and a range of categories, separates values
 * by category.  This function also guarantees that the order within each
 * category will be maintained.
 *
 * Note:  While the general convention of this library is to try to avoid
 * heap allocations whenever possible so that multithreaded code scales well and
 * false pointers aren't an issue, this function allocates like crazy
 * because there's basically no other way to implement it.  Don't use it in
 * performance-critical multithreaded code.
 *
 * Examples:
 * ---
 * uint[] values = [1,2,3,4,5,6,7,8];
 * bool[] categories = [false, false, false, false, true, true, true, true];
 * auto separated = byCategory(values, categories);
 * auto tResult = studentsTTest(separated.values);
 * ---
 */
ElementType!(V)[][ElementType!(C)] byCategory(V, C)(V values, C categories)
if(isInputRange!(V) && isInputRange!(C) && !is(ElementType!C == bool)) {
    alias ElementType!(V) EV;
    alias ElementType!(C) EC;

    Appender!(EV[])[EC] aa;
    while(!values.empty && !categories.empty) {
        scope(exit) {
            values.popFront();
            categories.popFront();
        }

        auto category = categories.front;
        auto ptr = category in aa;
        if(ptr is null) {
            aa[category] = typeof(aa[category]).init;
            ptr = category in aa;
        }

        ptr.put(values.front);
    }

    EV[][EC] ret;
    foreach(k, v; aa) {
        ret[k] = v.data;
    }

    return ret;
}

/**Special case implementation for when ElementType!C is boolean.*/
ElementType!(V)[][2] byCategory(V, C)(V values, C categories)
if(isInputRange!(V) && isInputRange!(C) && is(ElementType!C == bool)) {
    typeof(return) ret;

    static if(hasLength!V || hasLength!C) {
        // Preallocate.
        size_t zeroIndex, oneIndex;

        static if(!hasLength!V) {
            ret[0].length = categories.length;
            ret[1].length = values.length;
        } else static if(!hasLength!C) {
            ret[0].length = values.length;
            ret[1].length = values.length;
        } else {
            immutable len = min(categories.length, values.length);
            ret[0].length = len;
            ret[1].length = len;
        }

        while(!values.empty && !categories.empty) {
            scope(exit) {
                values.popFront();
                categories.popFront();
            }

            auto category = categories.front;
            if(category) {
                ret[1][oneIndex++] = values.front;
            } else {
                ret[0][zeroIndex++] = values.front;
            }
        }

        ret[0] = ret[0][0..zeroIndex];
        ret[1] = ret[1][0..oneIndex];
    } else {
        auto app0 = appender!(ElementType!(V)[])();
        auto app1 = appender!(ElementType!(V)[])();

        while(!values.empty && !categories.empty) {
            scope(exit) {
                values.popFront();
                categories.popFront();
            }

            auto category = categories.front;
            if(category) {
                app1.put(values.front);
            } else {
                app0.put(values.front);
            }
        }

        ret[0] = app0.data;
        ret[1] = app1.data;
    }

    return ret;
}

unittest {
    int[] nums = [1,2,3,4,5,6,7,8,9];
    int[] categories = [0,1,2,0,1,2,0,1,2];

    // The filter is just to prevent having a length.
    auto categories2 = filter!"a == a"(map!"a % 2 == 0"(nums));

    auto result = byCategory(nums, categories);
    assert(result[0] == [1,4,7]);
    assert(result[1] == [2,5,8]);
    assert(result[2] == [3,6,9]);

    auto res2 = byCategory(filter!"a == a"(nums), categories2);
    assert(res2[0] == [1,3,5,7,9]);
    assert(res2[1] == [2,4,6,8]);

    auto res3 = byCategory(nums,
        [false, true, false, true, false, true, false, true, false]);
    assert(res2 == res3);
}

/**Finds the area under the ROC curve (a curve with sensitivity on the Y-axis
 * and 1 - specificity on the X-axis).  This is a useful metric for
 * determining how well a test statistic discriminates between two classes.
 * The following assumptions are made in this implementation:
 *
 * 1.  For some cutoff value c and test statistic T, your decision rule is of
 *     the form "Class A if T > c, Class B if T < c".
 *
 * 2.  In the case of ties, i.e. if class A and class B both have an identical
 *     value, linear interpolation is used.  This is because changing the
 *     value of c infinitesimally will change both sensitivity and specificity
 *     in these cases.
 */
double auroc(R1, R2)(R1 classATs, R2 classBTs)
if(isNumeric!(ElementType!R1) && isNumeric!(ElementType!R2)) {
    auto alloc = newRegionAllocator();
    auto classA = alloc.array(classATs);
    auto classB = alloc.array(classBTs);
    qsort(classA);
    qsort(classB);

    // Start cutoff at -infinity, such that we get everything in class A, i.e.
    // perfect specificity, zero sensitivity.  We arbitrarily define class B
    // as our "positive" class.
    double tp = 0, tn = classA.length, fp = 0, fn = classB.length;
    double[2] lastPoint = 0;

    Unqual!(CommonType!(ElementType!R1, ElementType!R2)) currentVal;

    ElementType!R1 popA() {
        tn--;
        fp++;
        auto ret = classA.front;
        classA.popFront();
        return ret;
    }

    ElementType!R2 popB() {
        fn--;
        tp++;
        auto ret = classB.front;
        classB.popFront();
        return ret;
    }

    double area = 0;
    while(!classA.empty && !classB.empty) {
        if(classA.front < classB.front) {
            currentVal = popA();
        } else {
            currentVal = popB();
        }

        // Handle ties.
        while(!classA.empty && classA.front == currentVal) {
            popA();
        }

        while(!classB.empty && classB.front == currentVal) {
            popB();
        }

        double[2] curPoint;
        curPoint[0] = 1.0 - tn / (fp + tn);
        curPoint[1] = tp / (tp + fn);

        immutable xDist = curPoint[0] - lastPoint[0];
        area += xDist * lastPoint[1];  // Rectangular part.
        area += xDist * 0.5 * (curPoint[1] - lastPoint[1]);  // Triangular part.
        lastPoint[] = curPoint[];
    }

    if(classA.length > 0 && classB.length == 0) {
        // Then we already have a sensitivity of 1, move straight to the right
        // to the point (1, 1).

        immutable xDist = 1 - lastPoint[0];
        area += xDist * lastPoint[1];  // Rectangular part.
        area += xDist * 0.5 * (1 - lastPoint[1]);  // Triangular part.
    }

    return area;
}

unittest {
    // Values worked out by hand on paper.  If you don't believe me, work
    // them out yourself.
    assert(auroc([4,5,6], [1,2,3]) == 1);
    assert(approxEqual(auroc([8,6,7,5,3,0,9], [3,6,2,4,3,6]), 0.6904762));
    assert(approxEqual(auroc([2,7,1,8,2,8,1,8], [3,1,4,1,5,9,2,6]), 0.546875));
}

///
T sign(T)(T num) pure nothrow if(is(typeof(num < 0))) {
    if (num > 0) return 1;
    if (num < 0) return -1;
    return 0;
}

unittest {
    assert(sign(3.14159265)==1);
    assert(sign(-3)==-1);
    assert(sign(-2.7182818)==-1);
}

///
/*Values up to 9,999 are pre-calculated and stored in
 * an immutable global array, for performance.  After this point, the gamma
 * function is used, because caching would take up too much memory, and if
 * done lazily, would cause threading issues.*/
double logFactorial(ulong n) {
    //Input is uint, can't be less than 0, no need to check.
    if(n < staticFacTableLen) {
        return logFactorialTable[cast(size_t) n];
    } else return logGamma(n + 1);
}

unittest {
    // Cache branch.
    assert(cast(uint) round(exp(logFactorial(4)))==24);
    assert(cast(uint) round(exp(logFactorial(5)))==120);
    assert(cast(uint) round(exp(logFactorial(6)))==720);
    assert(cast(uint) round(exp(logFactorial(7)))==5040);
    assert(cast(uint) round(exp(logFactorial(3)))==6);
    // Gamma branch.
    assert(approxEqual(logFactorial(12000), 1.007175584216837e5, 1e-14));
    assert(approxEqual(logFactorial(14000), 1.196610688711534e5, 1e-14));
}

///Log of (n choose k).
double logNcomb(ulong n, ulong k)
in {
    assert(k <= n);
} body {
    if(n < k) return -double.infinity;
    //Extra parentheses increase numerical accuracy.
    return logFactorial(n) - (logFactorial(k) + logFactorial(n - k));
}

unittest {
    assert(cast(uint) round(exp(logNcomb(4,2)))==6);
    assert(cast(uint) round(exp(logNcomb(30,8)))==5852925);
    assert(cast(uint) round(exp(logNcomb(28,5)))==98280);
}

static if(size_t.sizeof == 4) {
    private enum MAX_PERM_LEN = 12;
} else {
    private enum MAX_PERM_LEN = 20;
}

/**
A struct that generates all possible permutations of a sequence.

Notes:

Permutations are output in undefined order.

The array returned by front is recycled across iterations.  To preserve
it across iterations, wrap this range using map!"a.dup" or
map!"a.idup".

Bugs:  Only supports iterating over up to size_t.max permutations.
This means the max permutation length is 12 on 32-bit machines, or 20
on 64-bit.  This was a conscious tradeoff to allow this range to have a
length of type size_t, since iterating over such huge permutation spaces
would be insanely slow anyhow.

Examples:
---
double[][] res;
auto perm = map!"a.dup"(Perm!(double)([1.0, 2.0, 3.0][]));
foreach(p; perm) {
     res ~= p;
}

auto sorted = sort(res);
assert(sorted.canFind([1.0, 2.0, 3.0]));
assert(sorted.canFind([1.0, 3.0, 2.0]));
assert(sorted.canFind([2.0, 1.0, 3.0]));
assert(sorted.canFind([2.0, 3.0, 1.0]));
assert(sorted.canFind([3.0, 1.0, 2.0]));
assert(sorted.canFind([3.0, 2.0, 1.0]));
assert(sorted.length == 6);
---
 */
struct Perm(T) {
private:

    // Optimization:  Since we know this thing can't get too big (there's
    // an dstatsEnforce statement for it in the c'tor), just use arrays of the max
    // possible size for stuff and store them inline, if it's all just bytes.
    static if(T.sizeof == 1) {
        T[MAX_PERM_LEN] perm;
    } else {
        T* perm;
    }

    // The length of this range.
    size_t nPerms;

    ubyte[MAX_PERM_LEN] Is;
    ubyte currentIndex;

    // The length of these arrays.  Stored once to minimize overhead.
    ubyte len;
    alias const(T)[] PermArray;

public:
    /**Generate permutations from an input range.
     * Create a duplicate of this sequence
     * so that the original sequence is not modified.*/
    this(U)(U input)
    if(isForwardRange!(U)) {

        static if(ElementType!(U).sizeof > 1) {
            auto arr = array(input);
            dstatsEnforce(arr.length <= MAX_PERM_LEN, text(
                "Can't iterate permutations of an array this long.  (Max length:  ",
                        MAX_PERM_LEN, ")"));
            len = cast(ubyte) arr.length;
            perm = arr.ptr;
        } else {
            foreach(elem; input) {
                dstatsEnforce(len < MAX_PERM_LEN, text(
                    "Can't iterate permutations of an array this long.  (Max length:  ",
                        MAX_PERM_LEN, ")"));

                perm[len++] = elem;
            }
        }

        popFront();

        nPerms = 1;
        for(size_t i = 2; i <= len; i++) {
            nPerms *= i;
        }
    }

    /**
    Returns the current permutation.  The array is const because it is
    recycled across iterations and modifying it would destroy the state of
    the permutation generator.
    */
    @property const(T)[] front() {
        return perm[0..len];
    }

    /**
    Get the next permutation in the sequence.  This will overwrite the
    contents of the array returned by the last call to front.
    */
    void popFront() {
        if(len == 0) {
            nPerms--;
            return;
        }
        if(currentIndex == len - 1) {
            currentIndex--;
            nPerms--;
            return;
        }

        uint max = len - currentIndex;
        if(Is[currentIndex] == max) {

            if(currentIndex == 0) {
                nPerms--;
                assert(nPerms == 0, to!string(nPerms));
                return;
            }

            Is[currentIndex..len] = 0;
            currentIndex--;
            return popFront();
        } else {
            rotateLeft(perm[currentIndex..len]);
            Is[currentIndex]++;
            currentIndex++;
            return popFront();
        }
    }

    ///
    @property bool empty() {
        return nPerms == 0;
    }

    /**
    The number of permutations left.
     */
    @property size_t length() const pure nothrow {
        return nPerms;
    }

    ///
    @property typeof(this) save() {
        auto ret = this;
        static if(T.sizeof > 1) {
            ret.perm = (ret.perm[0..len].dup).ptr;
        }
        return ret;
    }
}

/**
Create a Perm struct from a range or of a set of bounds.

Examples:
---
auto p = perm([1,2,3]);  // All permutations of [1,2,3].
auto p = perm(5);  // All permutations of [0,1,2,3,4].
auto p = perm(-1, 2); // All permutations of [-1, 0, 1].
---
 */
auto perm(T...)(T stuff) {
    static if(isForwardRange!(T[0])) {
        return Perm!(ElementType!T)(stuff);
    } else static if(T.length == 1) {
        static assert(isIntegral!(T[0]),
            "If one argument is passed to perm(), it must be an integer.");

        dstatsEnforce(stuff[0] >= 0, "Cannot generate permutations of length < 0.");
        dstatsEnforce(stuff[0] <= MAX_PERM_LEN, text(
            "Can't iterate permutations of an array of length ",
            stuff[0], ".  (Max length:  ", MAX_PERM_LEN, ")"));

        // Optimization:  Since we know the lower
        // bound is zero and the upper bound can't be > byte.max, use bytes
        // instead of bigger integer types.
        return Perm!byte(seq(cast(byte) 0, cast(byte) stuff[0]));
    } else {
        static assert(stuff.length == 2);
        return Perm!(CommonType!(T[0], T[1]))(seq(stuff[0], stuff[1]));
    }
}

unittest {
    // Test degenerate case of len == 0;
    uint nZero = 0;
    foreach(elem; perm(0)) {
        assert(elem.length == 0);
        nZero++;
    }
    assert(nZero == 1);

    double[][] res;
    auto p1 = map!"a.dup"(perm([1.0, 2.0, 3.0][]));
    assert(p1.length == 6);
    foreach(p; p1) {
        res ~= p;
    }
    auto sortedRes = sort(res);
    assert(sortedRes.contains([1.0, 2.0, 3.0]));
    assert(sortedRes.contains([1.0, 3.0, 2.0]));
    assert(sortedRes.contains([2.0, 1.0, 3.0]));
    assert(sortedRes.contains([2.0, 3.0, 1.0]));
    assert(sortedRes.contains([3.0, 1.0, 2.0]));
    assert(sortedRes.contains([3.0, 2.0, 1.0]));
    assert(res.length == 6);
    byte[][] res2;
    auto perm2 = map!"a.dup"(perm(3));
    foreach(p; perm2) {
        res2 ~= p;
    }
    auto sortedRes2 = sort(res2);
    assert(sortedRes2.contains( to!(byte[])([0, 1, 2])));
    assert(sortedRes2.contains( to!(byte[])([0, 2, 1])));
    assert(sortedRes2.contains( to!(byte[])([1, 0, 2])));
    assert(sortedRes2.contains( to!(byte[])([1, 2, 0])));
    assert(sortedRes2.contains( to!(byte[])([2, 0, 1])));
    assert(sortedRes2.contains( to!(byte[])([2, 1, 0])));
    assert(res2.length == 6);

    // Indirect tests:  If the elements returned are unique, there are N! of
    // them, and they contain what they're supposed to contain, the result is
    // correct.
    auto perm3 = perm(0U, 6U);
    bool[uint[]] table;
    foreach(p; perm3) {
        table[p.idup] = true;
    }
    assert(table.length == 720);
    foreach(elem, val; table) {
        assert(elem.dup.insertionSort == [0U, 1, 2, 3, 4, 5]);
    }
    auto perm4 = perm(5);
    bool[byte[]] table2;
    foreach(p; perm4) {
        table2[p.idup] = true;
    }
    assert(table2.length == 120);
    foreach(elem, val; table2) {
        assert(elem.dup.insertionSort == to!(byte[])([0, 1, 2, 3, 4]));
    }
}

/**
Generates every possible combination of r elements of the given sequence, or
array indices from zero to N, depending on which c'tor is called.  Uses
an input range interface.

Note:  The buffer that is returned by front is recycled across iterations.
To duplicate it instead, use map!"a.dup" or map!"a.idup".

Bugs:  Only supports iterating over up to size_t.max combinations.
This was a conscious tradeoff to allow this range to have a
length of type size_t, since iterating over such huge combination spaces
would be insanely slow anyhow.

Examples:
---
auto comb1 = map!"a.dup"(Comb!(uint)(5, 2));
uint[][] vals;
foreach(c; comb1) {
    vals ~= c;
}
auto sorted = sort(vals);
assert(sorted.canFind([0u,1]));
assert(sorted.canFind([0u,2]));
assert(sorted.canFind([0u,3]));
assert(sorted.canFind([0u,4]));
assert(sorted.canFind([1u,2]));
assert(sorted.canFind([1u,3]));
assert(sorted.canFind([1u,4]));
assert(sorted.canFind([2u,3]));
assert(sorted.canFind([2u,4]));
assert(sorted.canFind([3u,4]));
assert(sorted.length == 10);
---
 */
struct Comb(T) {
private:
    int N;
    int R;
    int diff;
    uint* pos;
    T* myArray;
    T* chosen;
    size_t _length;

    alias const(T)[] CombArray;

    void popFrontNum() {
        int index = R - 1;
        for(; index != -1 && pos[index] == diff + index; --index) {}
        if(index == -1) {
            _length--;
            return;
        }
        pos[index]++;
        for(uint i = index + 1; i < R; ++i) {
            pos[i] = pos[index] + i - index;
        }
        _length--;
    }

    void popFrontArray() {
        int index = R - 1;
        for(; index != -1 && pos[index] == diff + index; --index) {}
        if(index == -1) {
            _length--;
            return;
        }
        pos[index]++;
        chosen[index] = myArray[pos[index]];
        for(uint i = index + 1; i < R; ++i) {
            pos[i] = pos[index] + i - index;
            chosen[i] = myArray[pos[i]];
        }
        _length--;
    }

    void setLen() {
        // Used at construction.
        auto rLen = exp( logNcomb(N, R));
        dstatsEnforce(rLen < size_t.max, "Too many combinations.");
        _length = roundTo!size_t(rLen);
    }

public:

    /**
    Ctor to generate all possible combinations of array indices for a length r
    array.  This is a special-case optimization and is faster than simply
    using the other ctor to generate all length r combinations from
    seq(0, length).

    For efficiency, uint is used instead of size_t since, on a 64-bit system,
    generating all possible combinations of an array bigger than uint.max
    wouldn't be feasible anyhow.
    */
    static if(is(T == uint)) {
        this(uint n, uint r)
        in {
            assert(n >= r);
        } body {
            if(r > 0) {
                pos = (seq(0U, r)).ptr;
                pos[r - 1]--;
            }
            N = n;
            R = r;
            diff = N - R;
            popFront();
            setLen();
        }
    }

    /**General ctor.  array is a sequence from which to generate the
     * combinations.  r is the length of the combinations to be generated.*/
    this(T[] array, uint r) {
        if(r > 0) {
            pos = (seq(0U, r)).ptr;
            pos[r - 1]--;
        }
        N = to!uint(array.length);
        R = r;
        diff = N - R;
        auto temp = array.dup;
        myArray = temp.ptr;
        chosen = (new T[r]).ptr;
        foreach(i; 0..r) {
            chosen[i] = myArray[pos[i]];
        }
        popFront();
        setLen();
    }

    /**
    Gets the current combination.
    */
    @property const(T)[] front() {
        static if(!is(T == uint)) {
            return chosen[0..R].dup;
        } else {
            return (myArray is null) ? pos[0..R] : chosen[0..R];
        }
    }

    /**
    Advances to the next combination.  The array returned by front will be
    overwritten with the new results.
    */
    void popFront() {
        return (myArray is null) ? popFrontNum() : popFrontArray();
    }

    ///
    @property bool empty() const pure nothrow {
        return length == 0;
    }

    ///
    @property size_t length() const pure nothrow {
        return _length;
    }

    ///
    @property typeof(this) save() {
        auto ret = this;
        ret.pos = (pos[0..R].dup).ptr;

        if(chosen !is null) {
            ret.chosen = (chosen[0..R].dup).ptr;
        }

        return ret;
    }
}

/**Create a Comb struct from a range or of a set of bounds.
 *
 * Examples:
 * ---
 * auto c1 = comb([1,2,3], 2);  // Any two elements from [1,2,3].
 * auto c2 = comb(5, 3);  // Any three elements from [0,1,2,3,4].
 * ---
 */
auto comb(T)(T stuff, uint r) {
    static if(isForwardRange!(T)) {
        return Comb!(ElementType!T)(stuff, r);
    } else {
        static assert(isIntegral!T, "Can only call comb on ints and ranges.");
        return Comb!(uint)(cast(uint) stuff, r);
    }
}

unittest {
    // Test degenerate case of r == 0.  Shouldn't segfault.
    uint nZero = 0;
    foreach(elem; comb(5, 0)) {
        assert(elem.length == 0);
        nZero++;
    }
    assert(nZero == 1);

    nZero = 0;
    uint[] foo = [1,2,3,4,5];
    foreach(elem; comb(foo, 0)) {
        assert(elem.length == 0);
        nZero++;
    }
    assert(nZero == 1);

    // Test indexing verison first.
    auto comb1 = map!"a.dup"(comb(5, 2));
    uint[][] vals;
    foreach(c; comb1) {
        vals ~= c;
    }

    auto sortedVals = sort(vals);
    assert(sortedVals.contains([0u,1].dup));
    assert(sortedVals.contains([0u,2].dup));
    assert(sortedVals.contains([0u,3].dup));
    assert(sortedVals.contains([0u,4].dup));
    assert(sortedVals.contains([1u,2].dup));
    assert(sortedVals.contains([1u,3].dup));
    assert(sortedVals.contains([1u,4].dup));
    assert(sortedVals.contains([2u,3].dup));
    assert(sortedVals.contains([2u,4].dup));
    assert(sortedVals.contains([3u,4].dup));
    assert(vals.length == 10);

    // Now, test the array version.
    auto comb2 = map!"a.dup"(comb(seq(5U, 10U), 3));
    vals = null;
    foreach(c; comb2) {
        vals ~= c;
    }
    sortedVals = sort(vals);
    assert(sortedVals.contains([5u, 6, 7].dup));
    assert(sortedVals.contains([5u, 6, 8].dup));
    assert(sortedVals.contains([5u, 6, 9].dup));
    assert(sortedVals.contains([5u, 7, 8].dup));
    assert(sortedVals.contains([5u, 7, 9].dup));
    assert(sortedVals.contains([5u, 8, 9].dup));
    assert(sortedVals.contains([6U, 7, 8].dup));
    assert(sortedVals.contains([6u, 7, 9].dup));
    assert(sortedVals.contains([6u, 8, 9].dup));
    assert(sortedVals.contains([7u, 8, 9].dup));
    assert(vals.length == 10);

    // Now a test of a larger dataset where more subtle bugs could hide.
    // If the values returned are unique even after sorting, are composed of
    // the correct elements, and there is the right number of them, this thing
    // works.

    bool[uint[]] results;  // Keep track of how many UNIQUE items we have.
    auto comb3 = Comb!(uint)(seq(10U, 22U), 6);
    foreach(c; comb3) {
        auto dupped = c.dup.sort();
        // Make sure all elems are unique and within range.
        assert(dupped.length == 6);
        assert(dupped[0] > 9 && dupped[0] < 22);
        foreach(i; 1..dupped.length) {
            // Make sure elements are unique.  Remember, the array is sorted.
            assert(dupped[i] > dupped[i - 1]);
            assert(dupped[i] > 9 && dupped[i] < 22);
        }
        results[dupped.release().idup] = true;
    }
    assert(results.length == 924);  // (12 choose 6).
}

/**Converts a range with arbitrary element types (usually strings) to a
 * range of reals lazily.  Ignores any elements that could not be successfully
 * converted.  Useful for creating an input range that can be used with this
 * lib out of a File without having to read the whole file into an array first.
 * The advantages to this over just using std.algorithm.map are that it's
 * less typing and that it ignores non-convertible elements, such as blank
 * lines.
 *
 * If rangeIn is an inputRange, then the result of this function is an input
 * range.  Otherwise, the result is a forward range.
 *
 * Note:  The reason this struct doesn't have length or random access,
 * even if this is supported by rangeIn, is because it has to be able to
 * filter out non-convertible elements.
 *
 * Examples:
 * ---
 * // Perform a T-test to see whether the mean of the data being input as text
 * // from stdin is different from zero.  This data might not even fit in memory
 * // if it were read in eagerly.
 *
 * auto myRange = toNumericRange( stdin.byLine() );
 * TestRes result = studentsTTest(myRange);
 * writeln(result);
 * ---
 */
ToNumericRange!R toNumericRange(R)(R rangeIn) if(isInputRange!R) {
    alias ToNumericRange!R RT;
    return RT(rangeIn);
}

///
struct ToNumericRange(R) if(isInputRange!R) {
private:
    alias ElementType!R E;
    R inputRange;
    real _front;

public:
    this(R inputRange) {
        this.inputRange = inputRange;
        try {
            _front = to!real(inputRange.front);
        } catch(ConvException) {
            popFront();
        }
    }

    @property real front() {
        return _front;
    }

    void popFront() {
        while(true) {
            inputRange.popFront();
            if(inputRange.empty) {
                return;
            }
            auto inFront = inputRange.front;

            // If inFront is some string, strip the whitespace.
            static if( is(typeof(strip(inFront)))) {
                inFront = strip(inFront);
            }

            try {
                _front = to!real(inFront);
                return;
            } catch(ConvException) {
                continue;
            }
        }
    }

    @property bool empty() {
        return inputRange.empty;
    }
}

unittest {
    // Test both with non-convertible element as first element and without.
    // This is because non-convertible elements as the first element are
    // handled as a special case in the implementation.
    string[2] dataArr = ["3.14\n2.71\n8.67\nabracadabra\n362436",
                 "foobar\n3.14\n2.71\n8.67\nabracadabra\n362436"];

    foreach(data; dataArr) {
        std.file.write("NumericFileTestDeleteMe.txt", data);
        scope(exit) std.file.remove("NumericFileTestDeleteMe.txt");
        auto myFile = File("NumericFileTestDeleteMe.txt");
        auto rng = toNumericRange(myFile.byLine());
        assert(approxEqual(rng.front, 3.14));
        rng.popFront;
        assert(approxEqual(rng.front, 2.71));
        rng.popFront;
        assert(approxEqual(rng.front, 8.67));
        rng.popFront;
        assert(approxEqual(rng.front, 362435));
        assert(!rng.empty);
        rng.popFront;
        assert(rng.empty);
        myFile.close();
    }
}

// Used for ILP optimizations.
package template StaticIota(size_t upTo) {
    static if(upTo == 0) {
        alias TypeTuple!() StaticIota;
    } else {
        alias TypeTuple!(StaticIota!(upTo - 1), upTo - 1) StaticIota;
    }
}

package:

// If SciD is available, it is used for matrix storage and operations
// such as Cholesky decomposition.  Otherwise, some old, crappy routines
// that were around since before SciD are used.  
version(scid) {
    import scid.matvec, scid.linalg, scid.exception;
    
    alias ExternalMatrixView!double DoubleMatrix;
    
    DoubleMatrix doubleMatrix(
        size_t nRows, 
        size_t nColumns,
        RegionAllocator alloc
    ) {
        return DoubleMatrix(nRows, nColumns, alloc);
    }
    
    void choleskySolve(DoubleMatrix a, double[] b, double[] x) {
        choleskyDestructive(a);
        auto ans = externalVectorView(x);
        auto vec = externalVectorView(b);
        scid.linalg.choleskySolve(a, vec, ans);
    }
    
    void invert(DoubleMatrix from, DoubleMatrix to) {
        try {
            to[] = inv(from);
        } catch(Exception) {
            foreach(i; 0..to.rows) foreach(j; 0..to.columns) {
                to[i, j] = double.nan;
            }
        }
    }
        
} else {
    version = noscid;
    
    struct DoubleMatrix {
        double[][] arrayOfArrays;
        //alias arrayOfArraysFun this;
        
        const pure nothrow @property {
            size_t rows() {
                return arrayOfArrays.length;
            }
            
            size_t columns() {
                return (arrayOfArrays.length) ? 
                    (arrayOfArrays[0].length) : 0;
            }
        }                    
        
        ref double opIndex(size_t i, size_t j) {
            return arrayOfArrays[i][j];
        }
    }
    
    DoubleMatrix doubleMatrix(
        size_t nRows, 
        size_t nColumns,
        RegionAllocator alloc
    ) {
        return DoubleMatrix(
            alloc.uninitializedArray!(double[][])(nRows, nColumns)
        );
    }    
    
    void invert(DoubleMatrix from, DoubleMatrix to) {
        invert(from.arrayOfArrays, to.arrayOfArrays);
    }
    
    // Uses Gauss-Jordan elim. w/ row pivoting to invert from.  Stores the results
    // in to and leaves from in an undefined state.
    package void invert(double[][] from, double[][] to) {
        // Normalize.
        foreach(i, row; from) {
            double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..from.length]));
            row[] *= absMax;
            to[i][] = 0;
            to[i][i] = absMax;
        }

        foreach(col; 0..from.length) {
            size_t bestRow;
            double biggest = 0;
            foreach(row; col..from.length) {
                if(abs(from[row][col]) > biggest) {
                    bestRow = row;
                    biggest = abs(from[row][col]);
                }
            }

            swap(from[col], from[bestRow]);
            swap(to[col], to[bestRow]);
            immutable pivotFactor = from[col][col];

            foreach(row; 0..from.length) if(row != col) {
                immutable ratio = from[row][col] / pivotFactor;

                // If you're ever looking to optimize this code, good luck.  The
                // bottleneck is almost ENTIRELY this one line:
                from[row][] -= from[col][] * ratio;
                to[row][] -= to[col][] * ratio;
            }
        }

        foreach(i; 0..from.length) {
            immutable diagVal = from[i][i];
            from[i][] /= diagVal;
            to[i][] /= diagVal;
        }
    }

    unittest {
        auto mat = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]];
        auto toMat = [new double[3], new double[3], new double[3]];
        invert(mat, toMat);
        assert(approxEqual(toMat[0], [-0.625, -0.25, 0.375]));
        assert(approxEqual(toMat[1], [-0.25, 0.5, -0.25]));
        assert(approxEqual(toMat[2], [0.708333, -0.25, 0.041667]));
    }

    void solve(DoubleMatrix mat, double[] vec) {
        solve(mat.arrayOfArrays, vec);
    }

    // Solve a system of linear equations mat * b = vec for b.  The result is
    // stored in vec, and mat is left in an undefined state. Uses Gaussian
    // elimination w/ row pivoting.  Roughly 4x faster than using inversion to
    // solve the system, and uses roughly half the memory.
    void solve(double[][] mat, double[] vec) {
        // Normalize.
        foreach(i, row; mat) {
            double absMax = 1.0 / reduce!(max)(map!(abs)(row[0..mat.length]));
            row[] *= absMax;
            vec[i] *= absMax;
        }

        foreach(col; 0..mat.length - 1) {
            size_t bestRow;
            double biggest = 0;
            foreach(row; col..mat.length) {
                if(abs(mat[row][col]) > biggest) {
                    bestRow = row;
                    biggest = abs(mat[row][col]);
                }
            }

            swap(mat[col], mat[bestRow]);
            swap(vec[col], vec[bestRow]);
            immutable pivotFactor = mat[col][col];

            foreach(row; col + 1..mat.length) {
                immutable ratio = mat[row][col] / pivotFactor;

                // If you're ever looking to optimize this code, good luck.  The
                // bottleneck is almost ENTIRELY this one line:
                mat[row][col..$] -= mat[col][col..$] * ratio;
                vec[row] -= vec[col] * ratio;
            }
        }

        foreach(i; 0..mat.length) {
            double diagVal = mat[i][i];
            mat[i][] /= diagVal;
            vec[i] /= diagVal;
        }

        // Do back substitution.
        for(size_t row = mat.length - 1; row != size_t.max; row--) {
            auto v1 = vec[row + 1..$];
            auto v2 = mat[row][$ - v1.length..$];
            vec[row] -= dotProduct(v1, v2);
        }
    }

    unittest {
        auto mat = [[2.0, 1, -1], [-3.0, -1, 2], [-2.0, 1, 2]];
        auto vec = [8.0, -11, -3];
        solve(mat, vec);
        assert(approxEqual(vec, [2, 3, -1]));

        auto mat2 = [[1.0, 2, 3], [4.0, 7, 6], [7.0, 8, 9]];
        auto vec2 = [8.0, 6, 7];
        solve(mat2, vec2);
        assert(approxEqual(vec2, [-3.875, -0.75, 4.45833]));
    }

    // Cholesky decomposition functions adapted from Don Clugston's MathExtra
    // lib, used w/ permission.
    void choleskyDecompose(double[][] a, double[] diag) {
        immutable N = diag.length;

        foreach(i; 0..N) {
            const ai = a[i];
            double sum = ai[i];

            for(sizediff_t k = i - 1; k >= 0; --k) {
                sum -= ai[k] * ai[k];
            }

            if (sum > 0.0) {
                diag[i] = sqrt(sum);

                foreach(j; i + 1..N) {
                    sum = ai[j] - dotProduct(ai[0..i], a[j][0..i]);
                    a[j][i] = sum / diag[i];
                }
            } else {
                // not positive definite (could be caused by rounding errors)
                diag[i] = 0;
                // make this whole row zero so they have no further effect
                foreach(j; i + 1..N) a[j][i] = 0;
            }
        }
    }

    void choleskySolve(double[][] a, double[] diag, double[] b, double[] x) {
        immutable N = x.length;

        foreach(i; 0..N) {
            const ai = a[i];

            if(diag[i] > 0)  {
                double sum = b[i];
                sum -= dotProduct(ai[0..i], x[0..i]);
                x[i] = sum / diag[i];
            } else x[i] = 0; // skip pos definite rows
        }

        for(sizediff_t i = N - 1; i >= 0; --i) {
            if (diag[i] > 0) {
                double sum = x[i];
                for(sizediff_t k = i + 1; k < N; ++k) sum -= a[k][i] * x[k];
                x[i] = sum / diag[i];
            } else x[i] = 0; // skip pos definite rows
        }
        // Convert failed columns in solution to NANs if required.
        foreach(i; 0..N) {
            if(diag[i].isNaN() || diag[i] <= 0) x[i] = double.nan;
        }
    }

    void choleskySolve(DoubleMatrix a, double[] b, double[] x) {
        auto alloc = newRegionAllocator();
        auto diag = alloc.uninitializedArray!(double[])(x.length);
        choleskyDecompose(a.arrayOfArrays, diag);
        choleskySolve(a.arrayOfArrays, diag, b, x);
    }
}