/**A module for performing linear regression.  This module has an unusual
 * interface, as it is range-based instead of matrix based. Values for
 * independent variables are provided as either a tuple or a range of ranges.
 * This means that one can use, for example, map, to fit high order models and
 * lazily evaluate certain values.  (For details, see examples below.)
 *
 * 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.regress;

import std.math, std.traits, std.array, std.traits, std..string,
    std.exception, std.typetuple, std.typecons, std.numeric, std.parallelism;

import std.algorithm : map, copy, max, min, filter, reduce;

alias std.range.repeat repeat;

import dstats.alloc, std.range, std.conv, dstats.distrib, dstats.cor,
    dstats.base, dstats.summary, dstats.sort;

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

///
struct PowMap(ExpType, T)
if(isForwardRange!(T)) {
    Unqual!T range;
    Unqual!ExpType exponent;
    double cache;

    this(T range, ExpType exponent) {
        this.range = range;
        this.exponent = exponent;

        static if(isIntegral!ExpType) {
            cache = pow(cast(double) range.front, exponent);
        } else {
            cache = pow(cast(ExpType) range.front, exponent);
        }
    }

    @property double front() const pure nothrow {
        return cache;
    }

    void popFront() {
        range.popFront;
        if(!range.empty) {
            cache = pow(cast(double) range.front, exponent);
        }
    }

    @property typeof(this) save() {
        return this;
    }

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

/**Maps a forward range to a power determined at runtime.  ExpType is the type
 * of the exponent.  Using an int is faster than using a double, but obviously
 * less flexible.*/
PowMap!(ExpType, T) powMap(ExpType, T)(T range, ExpType exponent) {
    alias PowMap!(ExpType, T) RT;
    return RT(range, exponent);
}

// Very ad-hoc, does a bunch of matrix ops for linearRegress and
// linearRegressBeta.  Specifically, computes xTx and xTy.
// Written specifically to be efficient in the context used here.
private void rangeMatrixMulTrans(U, T...)(
    ref double[] xTy,
    ref DoubleMatrix xTx,
    U y, ref T matIn,
    RegionAllocator alloc
) {
    static if(isArray!(T[0]) &&
        isInputRange!(typeof(matIn[0][0])) && matIn.length == 1) {
        alias matIn[0] mat;
    } else {
        alias matIn mat;
    }

    bool someEmpty() {
        if(y.empty) {
            return true;
        }
        foreach(range; mat) {
            if(range.empty) {
                return true;
            }
        }
        return false;
    }

    void popAll() {
        foreach(ti, range; mat) {
            mat[ti].popFront;
        }
        y.popFront;
    }

    xTy[] = 0;
    xTx = doubleMatrix(mat.length, mat.length, alloc);
    foreach(i; 0..mat.length) foreach(j; 0..mat.length) {
        xTx[i, j] = 0;
    }

    auto alloc2 = newRegionAllocator();
    auto deltas = alloc2.uninitializedArray!(double[])(mat.length);

    // Using an algorithm similar to the one for Pearson cor to improve
    // numerical stability.  Calculate means and covariances, then
    // combine them:  Sum of squares = mean1 * N * mean2 + N * cov.
    double k = 0;
    double[] means = alloc2.uninitializedArray!(double[])(mat.length);
    means[] = 0;
    double yMean = 0;

    while(!someEmpty) {
        foreach(i, elem; mat) {
            deltas[i] = cast(double) elem.front;
        }

        immutable kMinus1 = k;
        immutable kNeg1 = 1 / ++k;
        deltas[] -= means[];
        means[] += deltas[] * kNeg1;

        immutable yDelta = cast(double) y.front - yMean;
        yMean += yDelta * kNeg1;

        foreach(i, delta1; deltas) {
            xTy[i] += kMinus1 * delta1 * kNeg1 * yDelta;
            xTx[i, i] += kMinus1 * delta1 * kNeg1 * delta1;

            foreach(j, delta2; deltas[0..i]) {
                xTx[i, j] += kMinus1 * delta1 * kNeg1 * delta2;
            }
        }
        popAll();
    }

    // k is n now that we're done looping over the data.
    alias k n;

    // mat now consists of covariance * n.    Divide by n and add mean1 * mean2
    // to get sum of products.
    foreach(i; 0..xTx.rows) foreach(j; 0..i + 1) {
        xTx[i, j] /= n;
        xTx[i, j] += means[i] * means[j];
    }

    // Similarly for the xTy vector
    foreach(i, ref elem; xTy) {
        xTy[i] /= n;
        elem += yMean * means[i];
    }
    symmetrize(xTx);
}

// Copies values from lower triangle to upper triangle.
private void symmetrize(ref DoubleMatrix mat) {
    foreach(i; 1..mat.rows) {
        foreach(j; 0..i) {
            mat[j, i] = mat[i, j];
        }
    }
}

/**Struct that holds the results of a linear regression.  It's a plain old
 * data struct.*/
struct RegressRes {
    /**The coefficients, one for each range in X.  These will be in the order
     * that the X ranges were passed in.*/
    double[] betas;

    /**The standard error terms of the X ranges passed in.*/
    double[] stdErr;

    /**The lower confidence bounds of the beta terms, at the confidence level
     * specificied.  (Default 0.95).*/
    double[] lowerBound;

    /**The upper confidence bounds of the beta terms, at the confidence level
     * specificied.  (Default 0.95).*/
    double[] upperBound;

    /**The P-value for the alternative that the corresponding beta value is
     * different from zero against the null that it is equal to zero.*/
    double[] p;

    /**The coefficient of determination.*/
    double R2;

    /**The adjusted coefficient of determination.*/
    double adjustedR2;

    /**The root mean square of the residuals.*/
    double residualError;

    /**The P-value for the model as a whole.  Based on an F-statistic.  The
     * null here is that the model has no predictive value, the alternative
     * is that it does.*/
    double overallP;

    // Just used internally.
    private static string arrStr(T)(T arr) {
        return text(arr)[1..$ - 1];
    }

    /**Print out the results in the default format.*/
    string toString() {
        return "Betas:  " ~ arrStr(betas) ~ "\nLower Conf. Int.:  " ~
            arrStr(lowerBound) ~ "\nUpper Conf. Int.:  " ~ arrStr(upperBound) ~
            "\nStd. Err:  " ~ arrStr(stdErr) ~ "\nP Values:  " ~ arrStr(p) ~
            "\nR^2:  " ~ text(R2) ~
            "\nAdjusted R^2:  " ~ text(adjustedR2) ~
            "\nStd. Residual Error:  " ~ text(residualError)
            ~ "\nOverall P:  " ~ text(overallP);
    }
}

/**Forward Range for holding the residuals from a regression analysis.*/
struct Residuals(F, U, T...) {
    static if(T.length == 1 && isForwardRange!(ElementType!(T[0]))) {
        alias T[0] R;
        alias typeof(array(R.init)) XType;
        enum bool needDup = true;
    } else {
        alias T R;
        alias staticMap!(Unqual, R) XType;
        enum bool needDup = false;
    }

    Unqual!U Y;
    XType X;
    F[] betas;
    double residual;
    bool _empty;

    void nextResidual() {
        double sum = 0;
        size_t i = 0;
        foreach(elem; X) {
            double frnt = elem.front;
            sum += frnt * betas[i];
            i++;
        }
        residual = Y.front - sum;
    }

    this(F[] betas, U Y, R X) {
        static if(is(typeof(X.length))) {
            dstatsEnforce(X.length == betas.length,
                "Betas and X must have same length for residuals.");
        } else {
            dstatsEnforce(walkLength(X) == betas.length,
                "Betas and X must have same length for residuals.");
        }

        static if(needDup) {
            this.X = array(X);
        } else {
            this.X = X;
        }

        foreach(i, elem; this.X) {
            static if(isForwardRange!(typeof(elem))) {
                this.X[i] = this.X[i].save;
            }
        }

        this.Y = Y;
        this.betas = betas;
        if(Y.empty) {
            _empty = true;
            return;
        }
        foreach(elem; X) {
            if(elem.empty) {
                _empty = true;
                return;
            }
        }
        nextResidual;
    }

    @property double front() const pure nothrow {
        return residual;
    }

    void popFront() {
        Y.popFront;
        if(Y.empty) {
            _empty = true;
            return;
        }
        foreach(ti, elem; X) {
            X[ti].popFront;
            if(X[ti].empty) {
                _empty = true;
                return;
            }
        }
        nextResidual;
    }

    @property bool empty() const pure nothrow {
        return _empty;
    }

    @property typeof(this) save() {
        auto ret = this;
        ret.Y = ret.Y.save;
        foreach(ti, xElem; ret.X) {
            ret.X[ti] = ret.X[ti].save;
        }

        return ret;
    }
}

/**Given the beta coefficients from a linear regression, and X and Y values,
 * returns a range that lazily computes the residuals.
 */
Residuals!(F, U, T) residuals(F, U, T...)(F[] betas, U Y, T X)
if(isFloatingPoint!F && isForwardRange!U && allSatisfy!(isForwardRange, T)) {
    alias Residuals!(F, U, T) RT;
    return RT(betas, Y, X);
}

// Compiles summary statistics while iterating, to allow ridge regression over
// input ranges.
private struct SummaryIter(R) {
    R range;
    MeanSD summ;

    this(R range) {
        this.range = range;
    }

    double front() @property {
        return range.front;
    }

    void popFront() {
        summ.put(range.front);
        range.popFront();
    }

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

    double mse() @property const pure nothrow { return summ.mse; }

    double N() @property const pure nothrow { return summ.N; }
}

private template SummaryType(R) {
    alias SummaryIter!R SummaryType;
}

/**
Perform a linear regression and return just the beta values.  The advantages
to just returning the beta values are that it's faster and that each range
needs to be iterated over only once, and thus can be just an input range.
The beta values are returned such that the smallest index corresponds to
the leftmost element of X.  X can be either a tuple or a range of input
ranges.  Y must be an input range.

If, after all X variables are passed in, a numeric type is passed as the last
parameter, this is treated as a ridge parameter and ridge regression is
performed.  Ridge regression is a form of regression that penalizes the L2 norm
of the beta vector and therefore results in more parsimonious models.
However, it makes statistical inference such as that supported by
linearRegress() difficult to impossible.  Therefore, linearRegress() doesn't
support ridges.

If no ridge parameter is passed, or equivalently if the ridge parameter is
zero, then ordinary least squares regression is performed.

Notes:  The X ranges are traversed in lockstep, but the traversal is stopped
at the end of the shortest one.  Therefore, using infinite ranges is safe.
For example, using repeat(1) to get an intercept term works.

References:

http://www.mathworks.com/help/toolbox/stats/ridge.html

Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S.
Fourth Edition. Springer, New York. ISBN 0-387-95457-0
(This is the citation for the MASS R package.)

Examples:
---
int[] nBeers = [8,6,7,5,3,0,9];
int[] nCoffees = [3,6,2,4,3,6,8];
int[] musicVolume = [3,1,4,1,5,9,2];
int[] programmingSkill = [2,7,1,8,2,8,1];
double[] betas = linearRegressBeta(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume));

// Now throw in a ridge parameter of 2.5.
double[] ridgeBetas = linearRegressBeta(programmingSkill, repeat(1), nBeers,
    nCoffees, musicVolume, map!"a * a"(musicVolume), 2.5);
---
 */
double[] linearRegressBeta(U, T...)(U Y, T XIn)
if(doubleInput!(U)) {
    double[] dummy;
    return linearRegressBetaBuf!(U, T)(dummy, Y, XIn);
}

/**
Same as linearRegressBeta, but allows the user to specify a buffer for
the beta terms.  If the buffer is too short, a new one is allocated.
Otherwise, the results are returned in the user-provided buffer.
 */
double[] linearRegressBetaBuf(U, TRidge...)(double[] buf, U Y, TRidge XRidge)
if(doubleInput!(U)) {
    auto alloc = newRegionAllocator();

    static if(isFloatingPoint!(TRidge[$ - 1]) || isIntegral!(TRidge[$ - 1])) {
        // ridge param.
        alias XRidge[$ - 1] ridgeParam;
        alias TRidge[0..$ - 1] T;
        alias XRidge[0..$ - 1] XIn;
        enum bool ridge = true;
        dstatsEnforce(ridgeParam >= 0,
            "Cannot do ridge regerssion with ridge param <= 0.");

        static SummaryIter!R summaryIter(R)(R range) {
            return typeof(return)(range);
        }

    } else {
        enum bool ridge = false;
        enum ridgeParam = 0;
        alias TRidge T;
        alias XRidge XIn;

        static R summaryIter(R)(R range) {
            return range;
        }
    }

    static if(isArray!(T[0]) && isInputRange!(typeof(XIn[0][0])) &&
        T.length == 1) {
        auto X = alloc.array(map!(summaryIter)(XIn[0]));
        alias typeof(X[0]) E;
    } else {
        static if(ridge) {
            alias staticMap!(SummaryType, T) XType;
            XType X;

            foreach(ti, elem; XIn) {
                X[ti] = summaryIter(elem);
            }
        } else {
            alias XIn X;
        }
    }

    DoubleMatrix xTx;
    double[] xTy = alloc.uninitializedArray!(double[])(X.length);
    double[] ret;
    if(buf.length < X.length) {
        ret = new double[X.length];
    } else {
        ret = buf[0..X.length];
    }

    rangeMatrixMulTrans(xTy, xTx, Y, X, alloc);

    static if(ridge) {
        if(ridgeParam > 0) {
            foreach(i, range; X) {
                // The whole matrix is scaled by N, as well as the xTy vector,
                // so scale the ridge param similarly.
                xTx[i, i] += ridgeParam * range.mse / range.N;
            }
        }
    }

    choleskySolve(xTx, xTy, ret);
    return ret;
}

/**
Perform a linear regression as in linearRegressBeta, but return a
RegressRes with useful stuff for statistical inference.  If the last element
of input is a real, this is used to specify the confidence intervals to
be calculated.  Otherwise, the default of 0.95 is used.  The rest of input
should be the elements of X.

When using this function, which provides several useful statistics useful
for inference, each range must be traversed twice.  This means:

1.  They have to be forward ranges, not input ranges.

2.  If you have a large amount of data and you're mapping it to some
    expensive function, you may want to do this eagerly instead of lazily.

Notes:

The X ranges are traversed in lockstep, but the traversal is stopped
at the end of the shortest one.  Therefore, using infinite ranges is safe.
For example, using repeat(1) to get an intercept term works.

If the confidence interval specified is exactly 0, this is treated as a
special case and confidence interval calculation is skipped.  This can speed
things up significantly and therefore can be useful in monte carlo and possibly
data mining contexts.

Bugs:  The statistical tests performed in this function assume that an
intercept term is included in your regression model.  If no intercept term
is included, the P-values, confidence intervals and adjusted R^2 values
calculated by this function will be wrong.

Examples:
---
int[] nBeers = [8,6,7,5,3,0,9];
int[] nCoffees = [3,6,2,4,3,6,8];
int[] musicVolume = [3,1,4,1,5,9,2];
int[] programmingSkill = [2,7,1,8,2,8,1];

// Using default confidence interval:
auto results = linearRegress(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume));

// Using user-specified confidence interval:
auto results = linearRegress(programmingSkill, repeat(1), nBeers, nCoffees,
    musicVolume, map!"a * a"(musicVolume), 0.8675309);
---
*/
RegressRes linearRegress(U, TC...)(U Y, TC input) {
    static if(is(TC[$ - 1] : double)) {
        double confLvl = input[$ - 1];
        enforceConfidence(confLvl);
        alias TC[0..$ - 1] T;
        alias input[0..$ - 1] XIn;
    } else {
        double confLvl = 0.95; // Default;
        alias TC T;
        alias input XIn;
    }

    auto alloc = newRegionAllocator();
    static if(isForwardRange!(T[0]) && isForwardRange!(typeof(XIn[0].front())) &&
        T.length == 1) {

        enum bool arrayX = true;
        alias typeof(XIn[0].front) E;
        E[] X = alloc.array(XIn[0]);
    } else static if(allSatisfy!(isForwardRange, T)) {
        enum bool arrayX = false;
        alias XIn X;
    } else {
        static assert(0, "Linear regression can only be performed with " ~
            "tuples of forward ranges or ranges of forward ranges.");
    }

    auto xTx = doubleMatrix(X.length, X.length, alloc);
    auto xTxNeg1 = doubleMatrix(X.length, X.length, alloc);
    double[] xTy = alloc.uninitializedArray!(double[])(X.length);

    static if(arrayX) {
        auto xSaved = alloc.array(X);
        foreach(ref elem; xSaved) {
            elem = elem.save;
        }
    } else {
        auto xSaved = X;
        foreach(ti, Type; X) {
            xSaved[ti] = X[ti].save;
        }
    }

    rangeMatrixMulTrans(xTy, xTx, Y.save, X, alloc);
    invert(xTx, xTxNeg1);
    auto betas = new double[X.length];
    foreach(i; 0..betas.length) {
        betas[i] = 0;
        foreach(j; 0..betas.length) {
            betas[i] += xTxNeg1[i, j] * xTy[j];
        }
    }

    X = xSaved;
    auto residuals = residuals(betas, Y, X);
    double S = 0;
    ulong n = 0;
    PearsonCor R2Calc;
    for(; !residuals.empty; residuals.popFront) {
        immutable residual = residuals.front;
        S += residual * residual;
        immutable Yfront = residuals.Y.front;
        immutable predicted = Yfront - residual;
        R2Calc.put(predicted, Yfront);
        n++;
    }
    immutable ulong df =  n - X.length;
    immutable double R2 = R2Calc.cor ^^ 2;
    immutable double adjustedR2 = 1.0L - (1.0L - R2) * ((n - 1.0L) / df);

    immutable double sigma2 = S / (n - X.length);

    double[] stdErr = new double[betas.length];
    foreach(i, ref elem; stdErr) {
        elem = sqrt( S * xTxNeg1[i, i] / df / n);
    }

    double[] lowerBound, upperBound;
    if(confLvl == 0) {
        // Then we're going to skip the computation to save time.  (See below.)
        lowerBound = betas;
        upperBound = betas;
    } else {
        lowerBound = new double[betas.length];
        upperBound = new double[betas.length];
    }
    auto p = new double[betas.length];

    foreach(i, beta; betas) {
        try {
            p[i] = 2 * min(studentsTCDF(beta / stdErr[i], df),
                           studentsTCDFR(beta / stdErr[i], df));
        } catch(DstatsArgumentException) {
            // Leave it as a NaN.
        }

        if(confLvl > 0) {
            // Skip confidence level computation if level is zero, to save
            // computation time.  This is important in monte carlo and possibly
            // data mining contexts.
            try {
                double delta = invStudentsTCDF(0.5 * (1 - confLvl), df) *
                     stdErr[i];
                upperBound[i] = beta - delta;
                lowerBound[i] = beta + delta;
            } catch(DstatsArgumentException) {
                // Leave confidence bounds as NaNs.
            }
        }
    }

    double F = (R2 / (X.length - 1)) / ((1 - R2) / (n - X.length));
    double overallP;
    try {
        overallP = fisherCDFR(F, X.length - 1, n - X.length);
    } catch(DstatsArgumentException) {
        // Leave it as a NaN.
    }

    return RegressRes(betas, stdErr, lowerBound, upperBound, p, R2,
        adjustedR2, sqrt(sigma2), overallP);
}


/**Struct returned by polyFit.*/
struct PolyFitRes(T) {

    /**The array of PowMap ranges created by polyFit.*/
    T X;

    /**The rest of the results.  This is alias this'd.*/
    RegressRes regressRes;
    alias regressRes this;
}

/**Convenience function that takes a forward range X and a forward range Y,
 * creates an array of PowMap structs for integer powers from 0 through N,
 * and calls linearRegressBeta.
 *
 * Returns:  An array of doubles.  The index of each element corresponds to
 * the exponent.  For example, the X<sup>2</sup> term will have an index of
 * 2.
 */
double[] polyFitBeta(T, U)(U Y, T X, uint N, double ridge = 0) {
    double[] dummy;
    return polyFitBetaBuf!(T, U)(dummy, Y, X, N);
}

/**Same as polyFitBeta, but allows the caller to provide an explicit buffer
 * to return the coefficients in.  If it's too short, a new one will be
 * allocated.  Otherwise, results will be returned in the user-provided buffer.
 */
double[] polyFitBetaBuf(T, U)(double[] buf, U Y, T X, uint N, double ridge = 0) {
    auto alloc = newRegionAllocator();
    auto pows = alloc.uninitializedArray!(PowMap!(uint, T)[])(N + 1);
    foreach(exponent; 0..N + 1) {
        pows[exponent] = powMap(X, exponent);
    }

    if(ridge == 0) {
        return linearRegressBetaBuf(buf, Y, pows);
    } else {
        return linearRegressBetaBuf(buf, Y, pows, ridge);
    }
}

/**Convenience function that takes a forward range X and a forward range Y,
 * creates an array of PowMap structs for integer powers 0 through N,
 * and calls linearRegress.
 *
 * Returns:  A PolyFitRes containing the array of PowMap structs created and
 * a RegressRes.  The PolyFitRes is alias this'd to the RegressRes.*/
PolyFitRes!(PowMap!(uint, T)[])
polyFit(T, U)(U Y, T X, uint N, double confInt = 0.95) {
    enforceConfidence(confInt);
    auto pows = new PowMap!(uint, T)[N + 1];
    foreach(exponent; 0..N + 1) {
        pows[exponent] = powMap(X, exponent);
    }
    alias PolyFitRes!(typeof(pows)) RT;
    RT ret;
    ret.X = pows;
    ret.regressRes = linearRegress(Y, pows, confInt);
    return ret;
}

unittest {
    // These are a bunch of values gleaned from various examples on the Web.
    double[] heights = [1.47,1.5,1.52,1.55,1.57,1.60,1.63,1.65,1.68,1.7,1.73,1.75,
        1.78,1.8,1.83];
    double[] weights = [52.21,53.12,54.48,55.84,57.2,58.57,59.93,61.29,63.11,64.47,
        66.28,68.1,69.92,72.19,74.46];
    float[] diseaseSev = [1.9, 3.1, 3.3, 4.8, 5.3, 6.1, 6.4, 7.6, 9.8, 12.4];
    ubyte[] temperature = [2,1,5,5,20,20,23,10,30,25];

    // Values from R.
    auto res1 = polyFit(diseaseSev, temperature, 1);
    assert(approxEqual(res1.betas[0], 2.6623));
    assert(approxEqual(res1.betas[1], 0.2417));
    assert(approxEqual(res1.stdErr[0], 1.1008));
    assert(approxEqual(res1.stdErr[1], 0.0635));
    assert(approxEqual(res1.p[0], 0.0419));
    assert(approxEqual(res1.p[1], 0.0052));
    assert(approxEqual(res1.R2, 0.644));
    assert(approxEqual(res1.adjustedR2, 0.6001));
    assert(approxEqual(res1.residualError, 2.03));
    assert(approxEqual(res1.overallP, 0.00518));


    auto res2 = polyFit(weights, heights, 2);
    assert(approxEqual(res2.betas[0], 128.813));
    assert(approxEqual(res2.betas[1], -143.162));
    assert(approxEqual(res2.betas[2], 61.960));

    assert(approxEqual(res2.stdErr[0], 16.308));
    assert(approxEqual(res2.stdErr[1], 19.833));
    assert(approxEqual(res2.stdErr[2], 6.008));

    assert(approxEqual(res2.p[0], 4.28e-6));
    assert(approxEqual(res2.p[1], 1.06e-5));
    assert(approxEqual(res2.p[2], 2.57e-7));

    assert(approxEqual(res2.R2, 0.9989, 0.0001));
    assert(approxEqual(res2.adjustedR2, 0.9987, 0.0001));

    assert(approxEqual(res2.lowerBound[0], 92.9, 0.01));
    assert(approxEqual(res2.lowerBound[1], -186.8, 0.01));
    assert(approxEqual(res2.lowerBound[2], 48.7, 0.01));
    assert(approxEqual(res2.upperBound[0], 164.7, 0.01));
    assert(approxEqual(res2.upperBound[1], -99.5, 0.01));
    assert(approxEqual(res2.upperBound[2], 75.2, 0.01));

    auto res3 = linearRegress(weights, repeat(1), heights, map!"a * a"(heights));
    assert(res2.betas == res3.betas);

    double[2] beta1Buf;
    auto beta1 = linearRegressBetaBuf
        (beta1Buf[], diseaseSev, repeat(1), temperature);
    assert(beta1Buf.ptr == beta1.ptr);
    assert(beta1Buf[] == beta1[]);
    assert(approxEqual(beta1, res1.betas));
    auto beta2 = polyFitBeta(weights, heights, 2);
    assert(approxEqual(beta2, res2.betas));

    auto res4 = linearRegress(weights, repeat(1), heights);
    assert(approxEqual(res4.p, 3.604e-14));
    assert(approxEqual(res4.betas, [-39.062, 61.272]));
    assert(approxEqual(res4.p, [6.05e-9, 3.60e-14]));
    assert(approxEqual(res4.R2, 0.9892));
    assert(approxEqual(res4.adjustedR2, 0.9884));
    assert(approxEqual(res4.residualError, 0.7591));
    assert(approxEqual(res4.lowerBound, [-45.40912, 57.43554]));
    assert(approxEqual(res4.upperBound, [-32.71479, 65.10883]));

    // Test residuals.
    assert(approxEqual(residuals(res4.betas, weights, repeat(1), heights),
        [1.20184170, 0.27367611,  0.40823237, -0.06993322,  0.06462305,
         -0.40354255, -0.88170814,  -0.74715188, -0.76531747, -0.63076120,
         -0.65892680, -0.06437053, -0.08253613,  0.96202014,  1.39385455]));

    // Test nonzero ridge parameters.
        // Values from R's MASS package.
    auto a = [1, 2, 3, 4, 5, 6, 7];
    auto b = [8, 6, 7, 5, 3, 0, 9];
    auto c = [2, 7, 1, 8, 2, 8, 1];

    // With a ridge param. of zero, ridge regression reduces to regular
    // OLS regression.
    assert(approxEqual(linearRegressBeta(a, repeat(1), b, c, 0),
        linearRegressBeta(a, repeat(1), b, c)));

    // Test the ridge regression. Values from R MASS package.
    auto ridge1 = linearRegressBeta(a, repeat(1), b, c, 1);
    auto ridge2 = linearRegressBeta(a, repeat(1), b, c, 2);
    auto ridge3 = linearRegressBeta(c, [[1,1,1,1,1,1,1], a, b], 10);
    assert(approxEqual(ridge1, [6.0357757, -0.2729671, -0.1337131]));
    assert(approxEqual(ridge2, [5.62367784, -0.22449854, -0.09775174]));
    assert(approxEqual(ridge3, [5.82653624, -0.05197246, -0.27185592 ]));
}

private MeanSD[] calculateSummaries(X...)(X xIn, RegionAllocator alloc) {
    // This is slightly wasteful because it sticks this shallow dup in
    // an unfreeable pos on TempAlloc.
    static if(X.length == 1 && isRoR!(X[0])) {
        auto ret = alloc.uninitializedArray!(MeanSD[])(xIn[0].length);
        auto alloc2 = newRegionAllocator();
        auto x = alloc2.array(xIn[0]);

        foreach(ref range; x) {
            range = range.save;
        }

        enum allHaveLength = hasLength!(ElementType!(typeof(x)));

    } else {
        auto ret = alloc2.uninitializedArray!MeanSD(xIn.length);
        alias xIn x;

        foreach(ti, R; X) {
            x[ti] = x[ti].save;
        }

        enum allHaveLength = allSatisfy!(hasLength, X);
    }

    ret[] = MeanSD.init;

    static if(allHaveLength) {
        size_t minLen = size_t.max;
        foreach(range; x) {
            minLen = min(minLen, range.length);
        }

        foreach(i, range; x) {
            ret[i] = meanStdev(take(range, minLen));
        }

    } else {
        bool someEmpty() {
            foreach(range; x) {
                if(range.empty) return true;
            }

            return false;
        }

        void popAll() {
            foreach(ti, elem; x) {
                x[ti].popFront();
            }
        }

        while(!someEmpty) {
            foreach(i, range; x) {
                ret[i].put(range.front);
            }
            popAll();
        }
    }

    return ret;
}

private double softThresh(double z, double gamma) {
    if(gamma >= abs(z)) {
        return 0;
    } else if(z > 0) {
        return z - gamma;
    } else {
        return z + gamma;
    }
}

private struct PreprocessedData {
    MeanSD[] xSumm;
    MeanSD ySumm;
    double[] y;
    double[][] x;
}

private PreprocessedData preprocessStandardize(Y, X...)
(Y yIn, X xIn, RegionAllocator alloc) {
    static if(X.length == 1 && isRoR!(X[0])) {
        auto xRaw = alloc.array(xIn[0]);
    } else {
        alias xIn xRaw;
    }

    auto summaries = calculateSummaries(xRaw, alloc);
    immutable minLen = to!size_t(
        reduce!min(
            map!"a.N"(summaries)
        )
    );

    auto x = alloc.uninitializedArray!(double[][])(summaries.length);
    foreach(i, range; xRaw) {
        x[i] = alloc.array(map!(to!double)(take(range, minLen)));
        x[i][] -= summaries[i].mean;
        x[i][] /= sqrt(summaries[i].mse);
    }

    double[] y;
    MeanSD ySumm;
    if(yIn.length) {
        y = alloc.array(map!(to!double)(take(yIn, minLen)));
        ySumm = meanStdev(y);
        y[] -= ySumm.mean;
    }

    return PreprocessedData(summaries, ySumm, y, x);
}

/**
Performs lasso (L1) and/or ridge (L2) penalized linear regression.  Due to the
way the data is standardized, no intercept term should be included in x
(unlike linearRegress and linearRegressBeta).  The intercept coefficient is
implicitly included and returned in the first element of the returned array.
Usage is otherwise identical.

Note:  Setting lasso equal to zero is equivalent to performing ridge regression.
       This can also be done with linearRegressBeta.  However, the
       linearRegressBeta algorithm is optimized for memory efficiency and
       large samples.  This algorithm is optimized for large feature sets.

Returns:  The beta coefficients for the regression model.

References:

Friedman J, et al Pathwise coordinate optimization. Ann. Appl. Stat.
2007;2:302-332.

Goeman, J. J., L1 penalized estimation in the Cox proportional hazards model.
Biometrical Journal 52(1), 70{84.

Eilers, P., Boer, J., Van Ommen, G., Van Houwelingen, H. 2001 Classification of
microarray data with penalized logistic regression. Proceedings of SPIE.
Progress in Biomedical Optics and Images vol. 4266, pp. 187-198

Douglas M. Hawkins and Xiangrong Yin.  A faster algorithm for ridge regression
of reduced rank data.  Computational Statistics & Data Analysis Volume 40,
Issue 2, 28 August 2002, Pages 253-262

http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
*/
double[] linearRegressPenalized(Y, X...)
(Y yIn, X xIn, double lasso, double ridge) {
    auto alloc = newRegionAllocator();

    auto preproc = preprocessStandardize(yIn, xIn, alloc);

    auto summaries = preproc.xSumm;
    auto ySumm = preproc.ySumm;
    auto x = preproc.x;
    auto y = preproc.y;

    auto betasFull = new double[x.length + 1];
    betasFull[] = 0;
    auto betas = betasFull[1..$];

    if(lasso > 0) {
        coordDescent(y, x, betas, lasso, ridge, null);
    } else if(y.length > x.length) {
        // Correct for different )#*$# scaling conventions.
        foreach(i, feature; x) {
            feature[] /= sqrt(summaries[i].mse);
        }

        linearRegressBetaBuf(betas, y, x, ridge);

        // More correction for diff. scaling factors.
        foreach(i, ref b; betas) {
            b /= sqrt(summaries[i].mse);
        }

    } else {
        ridgeLargeP(y, x, ridge, betas, null);
    }

    foreach(i, ref elem; betas) {
        elem /= sqrt(summaries[i].mse);
    }

    betasFull[0] = ySumm.mean;
    foreach(i, beta; betas) {
        betasFull[0] -= beta * summaries[i].mean;
    }

    return betasFull;
}

private void coordDescent(
    double[] y,
    double[][] x,
    double[] betas,
    double lasso,
    double ridge,
    double[] w
) {
    auto alloc = newRegionAllocator();

    auto predictions = alloc.uninitializedArray!(double[])(y.length);
    predictions[] = 0;

    void makePredictions() {
        foreach(j, beta; betas) {
            predictions[] += x[j][] * beta;
        }
    }

    if(reduce!max(0.0, map!abs(betas)) > 0) {
        makePredictions();
    }

    auto residuals = alloc.uninitializedArray!(double[])(y.length);

    uint iter = 0;
    enum maxIter = 10_000;
    enum relEpsilon = 1e-5;
    enum absEpsilon = 1e-10;
    immutable n = cast(double) y.length;
    auto perm = alloc.array(iota(0U, x.length));

    auto weightDots = alloc.uninitializedArray!(double[])(x.length);
    if(w.length == 0) {
        ridge /= n;
        lasso /= n;
        weightDots[] = 1;
    } else {
        foreach(j, col; x) {
            weightDots[j] = dotProduct(w, map!"a * a"(col));
        }
    }

    double doIter(double[] betas, double[][] x) {
        double maxRelError = 0;
        bool predictionsChanged = true;

        foreach(j, ref b; betas) {
            if(b == 0) {
                if(predictionsChanged) {
                    residuals[] = y[] - predictions[];
                }
            } else {
                predictions[] -= x[j][] * b;
                residuals[] = y[] - predictions[];
            }

            double z;
            if(w.length) {
                z = 0;
                foreach(i, weight; w) {
                    z += weight * x[j][i] * residuals[i];
                }
            } else {
                z = dotProduct(residuals, x[j]) / n;
            }

            immutable newB = softThresh(z, lasso) / (weightDots[j] + ridge);
            immutable absErr = abs(b - newB);
            immutable err = absErr / max(abs(b), abs(newB));

            if(absErr > absEpsilon) {
                maxRelError = max(maxRelError, err);
            }

            if(b == 0 && newB == 0) {
                predictionsChanged = false;
            } else {
                predictionsChanged = true;
            }

            b = newB;
            if(b != 0) {
                predictions[] += x[j][] * b;
            }
        }

        return maxRelError;
    }

    void toConvergence() {
        double maxRelErr = doIter(betas, x);
        iter++;
        if(maxRelErr < relEpsilon) return;

        static bool absGreater(double x, double y) { return abs(x) > abs(y); }

        while(iter < maxIter) {
            size_t split = 0;
            while(split < betas.length && abs(betas[split]) > 0) {
                split++;
            }

            try {
                qsort!absGreater(betas, x, perm, weightDots);
            } catch(SortException) {
                betas[] = double.nan;
                break;
            }

            maxRelErr = double.infinity;
            for(; !(maxRelErr < relEpsilon) && split < betas.length
            && iter < maxIter; iter++) {
                maxRelErr = doIter(betas[0..split], x[0..split]);
            }

            maxRelErr = doIter(betas, x);
            iter++;
            if(maxRelErr < relEpsilon) break;
        }
    }

    toConvergence();
    try {
        qsort(perm, x, betas);
    } catch(SortException) {
        return;
    }
}

/*
Ridge regression, case where P > N, where P is number of features and N
is number of samples.
*/
private void ridgeLargeP(
    const(double)[] y,
    const double[][] x,
    immutable double lambda,
    double[] betas,
    const double[] w = null
) {
    static if(haveSvd) {
        return eilersRidge(y, x, lambda, betas, w);
    } else {
        return shermanMorrisonRidge(y, x, lambda, betas, w);
    }
}

version(scid) {
    version(nodeps) {
        enum haveSvd = false;
    } else {
        enum haveSvd = true;
    }
} else {
    enum haveSvd = false;
}

/*
An implementation of ridge regression for large dimension.  Taken from:

Eilers, P., Boer, J., Van Ommen, G., Van Houwelingen, H. 2001 Classification of
microarray data with penalized logistic regression. Proceedings of SPIE.
Progress in Biomedical Optics and Images vol. 4266, pp. 187-198

This algorithm is very fast, O(N^2 * P) but requires singular value
decomposition.  Therefore, we only use if we're using SciD with full
dependency support.
*/
static if(haveSvd) private void eilersRidge(
    const(double)[] yArr,
    const double[][] x,
    immutable double lambda,
    double[] betas,
    const double[] weightArr
) {
    if(x.length == 0) return;
    auto alloc = newRegionAllocator();
    immutable n = x[0].length;
    immutable p = x.length;

    auto xMat = ExternalMatrixView!double(n, p, alloc);
    foreach(i; 0..n) foreach(j; 0..p) {
        xMat[i, j] = x[j][i];
    }

    typeof(svdDestructive(xMat, SvdType.economy, alloc)) svdRes;
    try {
        svdRes = svdDestructive(xMat, SvdType.economy, alloc);
    } catch(Exception) {
        betas[] = double.nan;
        return;
    }

    auto us = eval(svdRes.u * svdRes.s, alloc);
    ExternalMatrixView!double usWus;
    ExternalVectorView!double usWy;

    // Have to cast away const because ExternalVectorView doesn't play
    // nice w/ it yet.
    auto y = ExternalVectorView!double(cast(double[]) yArr);

    if(weightArr.length) {
        // Once we've multiplied s by u, we don't need it anymore.  Overwrite
        // its contents with the weight array to get weights in the form of
        // a diagonal matrix.
        auto w = svdRes.s;
        svdRes.s = typeof(svdRes.s).init;

        foreach(i, weight; weightArr) {
            w[i, i] = weight;
        }

        usWus = eval(us.t * w * us, alloc);
        usWy = eval(us.t * w * y, alloc);
    } else {
        usWus = eval(us.t * us, alloc);
        usWy = eval(us.t * y, alloc);
    }

    assert(usWus.rows == usWus.columns);
    assert(usWus.rows == n);
    foreach(i; 0..n) {
        usWus[i, i] += lambda;
    }

    auto theta = eval(inv(usWus) * usWy, alloc);

    // Work around weird SciD bug by transposing matrix manually.
    auto v = ExternalMatrixView!(double, StorageOrder.RowMajor)(
        svdRes.vt.columns,
        svdRes.vt.data[0..svdRes.vt.rows * svdRes.vt.columns]
    );
    eval(v * theta, betas);
}

/*
This algorithm is used as a fallback in case we don't have SVD available.
It's O(N P^2) instead of O(N^2 * P).  It was adapted from the
following paper:

Douglas M. Hawkins and Xiangrong Yin.  A faster algorithm for ridge regression
of reduced rank data.  Computational Statistics & Data Analysis Volume 40,
Issue 2, 28 August 2002, Pages 253-262

It also uses Wikipedia's page on Sherman-Morrison:

http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula

I simplified it by only doing rank-1 updates of inv(x' * x * lambda * I),
since I'm just trying to calculate the ridge estimate, not do efficient
cross-validation.  The idea is to start with an empty dataset.  Here,
inv(x' * x + lambda * I) = 1 / lambda * I.  Then, add each sample (row of x)
sequentially.  This is equivalent to adding s * s' to (x' * x + lambda * I)
where s is the vector of predictors for the sample.  This is known as a
dyadic product.  If p is the number of features/dimensions, and n is the
number of samples, we have n updates, each of which is O(P^2) for an
O(N * P^2) algorithm.  Naive algoriths would be O(P^3).
*/
static if(!haveSvd)
private void shermanMorrisonRidge(
    const(double)[] y,
    const(double[])[] x,  // Column major
    immutable double lambda,
    double[] betas,
    const(double)[] w
) in {
    assert(lambda > 0);
    foreach(col; x) assert(col.length == x[0].length);
    if(x.length) assert(y.length == x[0].length);
} body {
    auto alloc = newRegionAllocator();
    immutable p = x.length;
    if(p == 0) return;
    immutable n = x[0].length;

    auto v = alloc.uninitializedArray!(double[])(p);
    double[] u;
    if(w.length) {
        u = alloc.uninitializedArray!(double[])(p);
    } else {
        u = v;
    }

    auto vTxTxNeg1 = alloc.uninitializedArray!(double[])(p);
    auto xTxNeg1u = alloc.uninitializedArray!(double[])(p);

    // Before any updates are done, x'x = I * lambda, so inv(x'x) = I / lambda.
    auto xTxNeg1 = alloc.uninitializedArray!(double[][])(p, p);
    foreach(i; 0..p) foreach(j; 0..p) {
        xTxNeg1[i][j] = (i == j) ? (1.0 / lambda) : 0;
    }

    foreach(sampleIndex; 0..n) {
        copy(transversal(x[], sampleIndex), v);
        if(w.length) {
            u[] = w[sampleIndex] * v[];
        }

        // Calculate denominator of update:  1 + v' * xTxNeg1 * u
        vTxTxNeg1[] = 0;
        foreach(rowIndex, row; xTxNeg1) {
            vTxTxNeg1[] += v[rowIndex] * row[];
        }
        immutable denom = 1.0 + dotProduct(vTxTxNeg1, u);

        // Calculate numerator.  The parentheses indicate how the operation
        // is coded.  This is for computational efficiency.  Removing the
        // parentheses would be mathematically equivalent due to associativity:
        // (xTxNeg1 * u) * (v' * xTxNeg1).
        xTxNeg1u[] = 0;
        foreach(rowIndex, row; xTxNeg1) {
            xTxNeg1u[rowIndex] = dotProduct(row, u);
        }

        foreach(i, row; xTxNeg1) {
            immutable multiplier = xTxNeg1u[i] / denom;
            row[] -= multiplier * vTxTxNeg1[];
        }
    }

    auto xTy = alloc.uninitializedArray!(double[])(p);
    if(w.length) {
        auto yw = alloc.uninitializedArray!(double[])(n);
        yw[] = y[] * w[];
        y = yw;
    }

    foreach(colIndex, col; x) {
        xTy[colIndex] = dotProduct(col[], y[]);
    }

    foreach(rowIndex, row; xTxNeg1) {
        betas[rowIndex] = dotProduct(row, xTy);
    }
}

unittest {
    // Test ridge regression.  We have three impls for all kinds of diff.
    // scenarios.  See if they all agree.  Note that the ridiculously small but
    // nonzero lasso param is to force the use of the coord descent algo.
    auto y = new double[12];
    auto x = new double[][16];
    foreach(ref elem; x) elem = new double[12];
    x[0][] = 1;
    auto gen = Random(31415);  // For random but repeatable results.

    foreach(iter; 0..1000) {
        foreach(col; x[1..$]) foreach(ref elem; col) elem = rNormal(0, 1, gen);
        foreach(ref elem; y) elem = rNormal(0, 1, gen);
        immutable ridge = uniform(0.1, 10.0, gen);

        auto normalEq = linearRegressBeta(y, x, ridge);
        auto coordDescent = linearRegressPenalized(
            y, x[1..$], double.min_normal, ridge);
        auto linalgTrick = linearRegressPenalized(y, x[1..$], 0, ridge);

        // Every once in a blue moon coordinate descent doesn't converge that
        // well.  These small errors are of no practical significance, hence
        // the wide tolerance.  However, if the direct normal equations
        // and linalg trick don't agree extremely closely, then something's
        // fundamentally wrong.
        assert(approxEqual(normalEq, coordDescent, 0.02, 1e-4), text(
            normalEq, coordDescent));
        assert(approxEqual(linalgTrick, coordDescent, 0.02, 1e-4), text(
            linalgTrick, coordDescent));
        assert(approxEqual(normalEq, linalgTrick, 1e-6, 1e-8), text(
            normalEq, linalgTrick));
    }

    // Test stuff that's got some lasso in it.  Values from R's Penalized
    // package.
    y = [1.0, 2.0, 3, 4, 5, 6, 7];
    x = [[8.0, 6, 7, 5, 3, 0, 9],
         [3.0, 6, 2, 4, 3, 6, 8],
         [3.0, 1, 4, 1, 5, 9, 2],
         [2.0, 7, 1, 8, 2, 8, 1]];

    assert(approxEqual(linearRegressPenalized(y, x, 1, 0),
        [4.16316, -0.3603197, 0.6308278, 0, -0.2633263]));
    assert(approxEqual(linearRegressPenalized(y, x, 1, 3),
        [2.519590, -0.09116883, 0.38067757, 0.13122413, -0.05637939]));
    assert(approxEqual(linearRegressPenalized(y, x, 2, 0.1),
        [1.247235, 0, 0.4440735, 0.2023602, 0]));
    assert(approxEqual(linearRegressPenalized(y, x, 5, 7),
        [3.453787, 0, 0.10968736, 0.01253992, 0]));
}

/**
Computes a logistic regression using a maximum likelihood estimator
and returns the beta coefficients.  This is a generalized linear model with
the link function f(XB) = 1 / (1 + exp(XB)). This is generally used to model
the probability that a binary Y variable is 1 given a set of X variables.

For the purpose of this function, Y variables are interpreted as Booleans,
regardless of their type.  X may be either a range of ranges or a tuple of
ranges.  However, note that unlike in linearRegress, they are copied to an
array if they are not random access ranges.  Note that each value is accessed
several times, so if your range is a map to something expensive, you may
want to evaluate it eagerly.

If the last parameter passed in is a numeric value instead of a range,
it is interpreted as a ridge parameter and ridge regression is performed.  This
penalizes the L2 norm of the beta vector (in a scaled space) and results
in more parsimonious models.  It limits the usefulness of inference techniques
(p-values, confidence intervals), however, and is therefore not offered
in logisticRegres().

If no ridge parameter is passed, or equivalenty if the ridge parameter is
zero, then ordinary maximum likelihood regression is performed.

Note that, while this implementation of ridge regression was tested against
the R Design Package implementation, it uses slightly different conventions
that make the results not comparable without transformation.  dstats uses a
biased estimate of the variance to scale the beta vector penalties, while
Design uses an unbiased estimate.  Furthermore, Design penalizes by 1/2 of the
L2 norm, whereas dstats penalizes by the L2 norm.  Therefore, if n is the
sample size, and lambda is the penalty used with dstats, the proper penalty
to use in Design to get the same results is 2 * (n - 1) * lambda / n.

Also note that, as in linearRegress, repeat(1) can be used for the intercept
term.

Returns:  The beta coefficients for the regression model.

References:

http://en.wikipedia.org/wiki/Logistic_regression

http://socserv.mcmaster.ca/jfox/Courses/UCLA/logistic-regression-notes.pdf

S. Le Cessie and J. C. Van Houwelingen.  Ridge Estimators in Logistic
Regression.  Journal of the Royal Statistical Society. Series C
(Applied Statistics), Vol. 41, No. 1(1992), pp. 191-201

Frank E Harrell Jr (2009). Design: Design Package. R package version 2.3-0.
http://CRAN.R-project.org/package=Design
 */
double[] logisticRegressBeta(T, U...)(T yIn, U xRidge) {
    static if(isFloatingPoint!(U[$ - 1]) || isIntegral!(U[$ - 1])) {
        alias xRidge[$ - 1] ridge;
        alias xRidge[0..$ - 1] xIn;
    } else {
        enum double ridge = 0.0;
        alias xRidge xIn;
    }

    return logisticRegressImpl(false, ridge, yIn, xIn).betas;
}

/**
Plain old data struct to hold the results of a logistic regression.
*/
struct LogisticRes {
    /**The coefficients, one for each range in X.  These will be in the order
     * that the X ranges were passed in.*/
    double[] betas;

    /**The standard error terms of the X ranges passed in.*/
    double[] stdErr;

    /**
    The Wald lower confidence bounds of the beta terms, at the confidence level
    specificied.  (Default 0.95).*/
    double[] lowerBound;

    /**
    The Wald upper confidence bounds of the beta terms, at the confidence level
    specificied.  (Default 0.95).*/
    double[] upperBound;

    /**
    The P-value for the alternative that the corresponding beta value is
    different from zero against the null that it is equal to zero.  These
    are calculated using the Wald Test.*/
    double[] p;

    /**
    The log likelihood for the null model.
    */
    double nullLogLikelihood;

    /**
    The log likelihood for the model fit.
    */
    double logLikelihood;

    /**
    Akaike Information Criterion, which is a complexity-penalized goodness-
    of-fit score, equal to 2 * k - 2 log(L) where L is the log likelihood and
    k is the number of parameters.
    */
    double aic() const pure nothrow @property @safe {
        return 2 * (betas.length - logLikelihood);
    }

    /**
    The P-value for the model as a whole, based on the likelihood ratio test.
    The null here is that the model has no predictive value, the alternative
    is that it does have predictive value.*/
    double overallP;

    // Just used internally.
    private static string arrStr(T)(T arr) {
        return text(arr)[1..$ - 1];
    }

    /**Print out the results in the default format.*/
    string toString() {
        return "Betas:  " ~ arrStr(betas) ~ "\nLower Conf. Int.:  " ~
            arrStr(lowerBound) ~ "\nUpper Conf. Int.:  " ~ arrStr(upperBound) ~
            "\nStd. Err:  " ~ arrStr(stdErr) ~ "\nP Values:  " ~ arrStr(p) ~
            "\nNull Log Likelihood:  " ~ text(nullLogLikelihood) ~
            "\nLog Likelihood:  " ~ text(logLikelihood) ~
            "\nAIC:  " ~ text(aic) ~
            "\nOverall P:  " ~ text(overallP);
    }
}

/**
Similar to logisticRegressBeta, but returns a LogisticRes with useful stuff for
statistical inference.  If the last element of input is a floating point
number instead of a range, it is used to specify the confidence interval
calculated.  Otherwise, the default of 0.95 is used.

References:

http://en.wikipedia.org/wiki/Wald_test
http://en.wikipedia.org/wiki/Akaike_information_criterion
*/
LogisticRes logisticRegress(T, V...)(T yIn, V input) {
    return logisticRegressImpl!(T, V)(true, 0, yIn, input);
}

unittest {
    // Values from R.  Confidence intervals from confint.default().
    // R doesn't automatically calculate likelihood ratio P-value, and reports
    // deviations instead of log likelihoods.  Deviations are just
    // -2 * likelihood.
    alias approxEqual ae;  // Save typing.

    // Start with the basics, with X as a ror.
    auto y1 =  [1,   0, 0, 0, 1, 0, 0];
    auto x1 = [[1.0, 1 ,1 ,1 ,1 ,1 ,1],
              [8.0, 6, 7, 5, 3, 0, 9]];
    auto res1 = logisticRegress(y1, x1);
    assert(ae(res1.betas[0], -0.98273));
    assert(ae(res1.betas[1], 0.01219));
    assert(ae(res1.stdErr[0], 1.80803));
    assert(ae(res1.stdErr[1], 0.29291));
    assert(ae(res1.p[0], 0.587));
    assert(ae(res1.p[1], 0.967));
    assert(ae(res1.aic, 12.374));
    assert(ae(res1.logLikelihood, -0.5 * 8.3758));
    assert(ae(res1.nullLogLikelihood, -0.5 * 8.3740));
    assert(ae(res1.lowerBound[0], -4.5264052));
    assert(ae(res1.lowerBound[1], -0.5618933));
    assert(ae(res1.upperBound[0], 2.560939));
    assert(ae(res1.upperBound[1], 0.586275));

    // Use tuple.
    auto y2   = [1,0,1,1,0,1,0,0,0,1,0,1];
    auto x2_1 = [3,1,4,1,5,9,2,6,5,3,5,8];
    auto x2_2 = [2,7,1,8,2,8,1,8,2,8,4,5];
    auto res2 = logisticRegress(y2, repeat(1), x2_1, x2_2);
    assert(ae(res2.betas[0], -1.1875));
    assert(ae(res2.betas[1], 0.1021));
    assert(ae(res2.betas[2], 0.1603));
    assert(ae(res2.stdErr[0], 1.5430));
    assert(ae(res2.stdErr[1], 0.2507));
    assert(ae(res2.stdErr[2], 0.2108));
    assert(ae(res2.p[0], 0.442));
    assert(ae(res2.p[1], 0.684));
    assert(ae(res2.p[2], 0.447));
    assert(ae(res2.aic, 21.81));
    assert(ae(res2.nullLogLikelihood, -0.5 * 16.636));
    assert(ae(res2.logLikelihood, -0.5 * 15.810));
    assert(ae(res2.lowerBound[0], -4.2116584));
    assert(ae(res2.lowerBound[1], -0.3892603));
    assert(ae(res2.lowerBound[2], -0.2528110));
    assert(ae(res2.upperBound[0], 1.8366823));
    assert(ae(res2.upperBound[1], 0.5934631));
    assert(ae(res2.upperBound[2], 0.5733693));

    auto x2Intercept = [1,1,1,1,1,1,1,1,1,1,1,1];
    auto res2a = logisticRegress(y2,
        filter!"a.length"([x2Intercept, x2_1, x2_2]));
    foreach(ti, elem; res2a.tupleof) {
        assert(ae(elem, res2.tupleof[ti]));
    }

    // Use a huge range of values to test numerical stability.

    // The filter is to make y3 a non-random access range.
    auto y3 = filter!"a < 2"([1,1,1,1,0,0,0,0]);
    auto x3_1 = filter!"a > 0"([1, 1e10, 2, 2e10, 3, 3e15, 4, 4e7]);
    auto x3_2 = [1e8, 1e6, 1e7, 1e5, 1e3, 1e0, 1e9, 1e11];
    auto x3_3 = [-5e12, 5e2, 6e5, 4e3, -999999, -666, -3e10, -2e10];
    auto res3 = logisticRegress(y3, repeat(1), x3_1, x3_2, x3_3, 0.99);
    assert(ae(res3.betas[0], 1.115e0));
    assert(ae(res3.betas[1], -4.674e-15));
    assert(ae(res3.betas[2], -7.026e-9));
    assert(ae(res3.betas[3], -2.109e-12));
    assert(ae(res3.stdErr[0], 1.158));
    assert(ae(res3.stdErr[1], 2.098e-13));
    assert(ae(res3.stdErr[2], 1.878e-8));
    assert(ae(res3.stdErr[3], 4.789e-11));
    assert(ae(res3.p[0], 0.336));
    assert(ae(res3.p[1], 0.982));
    assert(ae(res3.p[2], 0.708));
    assert(ae(res3.p[3], 0.965));
    assert(ae(res3.aic, 12.544));
    assert(ae(res3.nullLogLikelihood, -0.5 * 11.0904));
    assert(ae(res3.logLikelihood, -0.5 * 4.5442));
    // Not testing confidence intervals b/c they'd be so buried in numerical
    // fuzz.


    // Test with a just plain huge dataset that R chokes for several minutes
    // on.  If you think this unittest is slow, try getting the reference
    // values from R.
    auto y4 = chain(
                take(cycle([0,0,0,0,1]), 500_000),
                take(cycle([1,1,1,1,0]), 500_000));
    auto x4_1 = iota(0, 1_000_000);
    auto x4_2 = array(map!(compose!(exp, "a / 1_000_000.0"))(x4_1));
    auto x4_3 = take(cycle([1,2,3,4,5]), 1_000_000);
    auto x4_4 = take(cycle([8,6,7,5,3,0,9]), 1_000_000);
    auto res4 = logisticRegress(y4, repeat(1), x4_1, x4_2, x4_3, x4_4, 0.99);
    assert(ae(res4.betas[0], -1.574));
    assert(ae(res4.betas[1], 5.625e-6));
    assert(ae(res4.betas[2], -7.282e-1));
    assert(ae(res4.betas[3], -4.381e-6));
    assert(ae(res4.betas[4], -8.343e-6));
    assert(ae(res4.stdErr[0], 3.693e-2));
    assert(ae(res4.stdErr[1], 7.188e-8));
    assert(ae(res4.stdErr[2], 4.208e-2));
    assert(ae(res4.stdErr[3], 1.658e-3));
    assert(ae(res4.stdErr[4], 8.164e-4));
    assert(ae(res4.p[0], 0));
    assert(ae(res4.p[1], 0));
    assert(ae(res4.p[2], 0));
    assert(ae(res4.p[3], 0.998));
    assert(ae(res4.p[4], 0.992));
    assert(ae(res4.aic, 1089339));
    assert(ae(res4.nullLogLikelihood, -0.5 * 1386294));
    assert(ae(res4.logLikelihood, -0.5 * 1089329));
    assert(ae(res4.lowerBound[0], -1.668899));
    assert(ae(res4.lowerBound[1], 5.439787e-6));
    assert(ae(res4.lowerBound[2], -0.8366273));
    assert(ae(res4.lowerBound[3], -4.27406e-3));
    assert(ae(res4.lowerBound[4], -2.111240e-3));
    assert(ae(res4.upperBound[0], -1.478623));
    assert(ae(res4.upperBound[1], 5.810089e-6));
    assert(ae(res4.upperBound[2], -6.198418e-1));
    assert(ae(res4.upperBound[3], 4.265302e-3));
    assert(ae(res4.upperBound[4], 2.084554e-3));

    // Test ridge stuff.
    auto ridge2 = logisticRegressBeta(y2, repeat(1), x2_1, x2_2, 3);
    assert(ae(ridge2[0], -0.40279319));
    assert(ae(ridge2[1], 0.03575638));
    assert(ae(ridge2[2], 0.05313875));

    auto ridge2_2 = logisticRegressBeta(y2, repeat(1), x2_1, x2_2, 2);
    assert(ae(ridge2_2[0], -0.51411490));
    assert(ae(ridge2_2[1], 0.04536590));
    assert(ae(ridge2_2[2], 0.06809964));
}

/// The logistic function used in logistic regression.
double logistic(double xb) pure nothrow @safe {
    return 1.0 / (1 + exp(-xb));
}

/**
Performs lasso (L1) and/or ridge (L2) penalized logistic regression.  Due to the
way the data is standardized, no intercept term should be included in x
(unlike logisticRegress and logisticRegressBeta).  The intercept coefficient is
implicitly included and returned in the first element of the returned array.
Usage is otherwise identical.

Note:  Setting lasso equal to zero is equivalent to performing ridge regression.
       This can also be done with logisticRegressBeta.  However, the
       logisticRegressBeta algorithm is optimized for memory efficiency and
       large samples.  This algorithm is optimized for large feature sets.

Returns:  The beta coefficients for the regression model.

References:

Friedman J, et al Pathwise coordinate optimization. Ann. Appl. Stat.
2007;2:302-332.

Goeman, J. J., L1 penalized estimation in the Cox proportional hazards model.
Biometrical Journal 52(1), 70{84.

Eilers, P., Boer, J., Van Ommen, G., Van Houwelingen, H. 2001 Classification of
microarray data with penalized logistic regression. Proceedings of SPIE.
Progress in Biomedical Optics and Images vol. 4266, pp. 187-198

Douglas M. Hawkins and Xiangrong Yin.  A faster algorithm for ridge regression
of reduced rank data.  Computational Statistics & Data Analysis Volume 40,
Issue 2, 28 August 2002, Pages 253-262

http://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
*/
double[] logisticRegressPenalized(Y, X...)
(Y yIn, X xIn, double lasso, double ridge) {
    auto alloc = newRegionAllocator();

    static assert(!isInfinite!Y, "Can't do regression with infinite # of Y's.");
    static if(isRandomAccessRange!Y) {
        alias yIn y;
    } else {
        auto y = toBools(yIn);
    }

    static if(X.length == 1 && isRoR!X) {
        enum bool tupleMode = false;
        static if(isForwardRange!X) {
            auto x = toRandomAccessRoR(y.length, xIn, alloc);
        } else {
            auto x = toRandomAccessRoR(y.length, alloc.array(xIn), alloc);
        }
    } else {
        enum bool tupleMode = true;
        auto xx = toRandomAccessTuple(xIn, alloc);
        auto x = xx.expand;
    }

    auto betas = new double[x.length + 1];
    if(y.length >= x.length && lasso == 0) {
        // Add intercept term.
        static if(tupleMode) {
            doMLENewton(betas, (double[]).init, ridge, y, repeat(1), x);
        } else static if(is(x == double[][])) {
            auto xInt = alloc.uninitializedArray!(double[])(betas.length);
            xInt[1..$] = x[];
            xInt[0] = alloc.array(repeat(1.0, y.length));
            doMLENewton(betas, (double[]).init, ridge, y, xInt);
        } else {
            // No choice but to dup the whole thing.
            auto xInt = alloc.uninitializedArray!(double[][])(betas.length);
            xInt[0] = alloc.array(repeat(1.0, y.length));

            foreach(i, ref arr; xInt[1..$]) {
                arr = alloc.array(
                    map!(to!double)(take(x[i], y.length))
                );
            }

            doMLENewton(betas, (double[]).init, ridge, y, xInt);
        }

    } else {
        logisticRegressPenalizedImpl(betas, lasso, ridge, y, x);
    }

    return betas;
}

unittest {
    // Test ridge regression.  We have three impls for all kinds of diff.
    // scenarios.  See if they all agree.  Note that the ridiculously small but
    // nonzero lasso param is to force the use of the coord descent algo.
    auto y = new bool[12];
    auto x = new double[][16];
    foreach(ref elem; x) elem = new double[12];
    x[0][] = 1;
    auto gen = Random(31415);  // For random but repeatable results.

    foreach(iter; 0..1000) {
        foreach(col; x[1..$]) foreach(ref elem; col) elem = rNormal(0, 1, gen);

        // Nothing will converge if y is all true's or all false's.
        size_t trueCount;
        do {
            foreach(ref elem; y) elem = cast(bool) rBernoulli(0.5, gen);
            trueCount = count!"a"(y);
        } while(trueCount == 0 || trueCount == y.length);

        immutable ridge = uniform(0.1, 10.0, gen);

        auto normalEq = logisticRegressBeta(y, x, ridge);
        auto coordDescent = logisticRegressPenalized(
            y, x[1..$], double.min_normal, ridge);
        auto linalgTrick = logisticRegressPenalized(y, x[1..$], 0, ridge);

        // Every once in a blue moon coordinate descent doesn't converge that
        // well.  These small errors are of no practical significance, hence
        // the wide tolerance.  However, if the direct normal equations
        // and linalg trick don't agree extremely closely, then something's
        // fundamentally wrong.
        assert(approxEqual(normalEq, coordDescent, 0.02, 1e-4), text(
            normalEq, coordDescent));
        assert(approxEqual(linalgTrick, coordDescent, 0.02, 1e-4), text(
            linalgTrick, coordDescent));
        assert(approxEqual(normalEq, linalgTrick, 1e-6, 1e-8), text(
            normalEq, linalgTrick));
    }

    assert(approxEqual(logisticRegressBeta(y, x[0], x[1], x[2]),
        logisticRegressPenalized(y, x[1], x[2], 0, 0)));
    assert(approxEqual(logisticRegressBeta(y, [x[0], x[1], x[2]]),
        logisticRegressPenalized(y, [x[1], x[2]], 0, 0)));
    assert(approxEqual(logisticRegressBeta(y, [x[0], x[1], x[2]]),
        logisticRegressPenalized(y,
        [to!(float[])(x[1]), to!(float[])(x[2])], 0, 0)));

    // Make sure the adding intercept stuff is right for the Newton path.
    //assert(logisticRegressBeta(x[0], x[1], x[2]) ==

    // Test stuff that's got some lasso in it.  Values from R's Penalized
    // package.
    y = [1, 0, 0, 1, 1, 1, 0];
    x = [[8.0, 6, 7, 5, 3, 0, 9],
         [3.0, 6, 2, 4, 3, 6, 8],
         [3.0, 1, 4, 1, 5, 9, 2],
         [2.0, 7, 1, 8, 2, 8, 1]];

    // Values from R's Penalized package.  Note that it uses a convention for
    // the ridge parameter such that Penalized ridge = 2 * dstats ridge.
    assert(approxEqual(logisticRegressPenalized(y, x, 1, 0),
        [1.642080, -0.22086515, -0.02587546,  0.00000000, 0.00000000 ]));
    assert(approxEqual(logisticRegressPenalized(y, x, 1, 3),
        [0.5153373, -0.04278257, -0.00888014,  0.01316831,  0.00000000]));
    assert(approxEqual(logisticRegressPenalized(y, x, 2, 0.1),
        [0.2876821, 0, 0., 0., 0]));
    assert(approxEqual(logisticRegressPenalized(y, x, 1.2, 7),
        [0.367613 , -0.017227631, 0.000000000, 0.003875104, 0.000000000]));
}

/**
This function performs loess regression.  Loess regression is a local
regression procedure, where a prediction of the dependent (y) variable
is made from an observation of the independent (x) variable by weighted
least squares over x values in the neighborhood of the value being evaluated.

In the future a separate function may be included to perform loess regression
with multiple predictors.  However, one predictor is the much more common
case and the multiple predictor case will require a much different API
and implementation, so for now only one predictor is supported.

Params:

y      = Observations of the dependent variable.
x      = Observations of the independent variable.
span   = The fraction of x observations considered to be "in the neighborhood"
         when performing local regression to predict the y value for a new x.
         For example, if 8 observations are provided and span == 0.5,
         the 4 nearest neighbors will be used on evaluation.
degree = The polynomial degree of the local regression.  Must be less than
         the number of neighbors (span * x.length).

Returns:

A Loess1D object.  This object can be used to make predictions based on
the loess model.  The computations involved are done lazily, i.e. this
function sets up the Loess1D instance and returns without computing
any regression models.

Examples:
---
auto x = [1, 2, 3, 4, 5, 6, 7];
auto y = [3, 6, 2, 4, 3, 6, 8];

// Build the Loess1D object.
auto model = loess1D(y, x, 0.5);

// Compute the weights for robust regression.
model.computeRobustWeights(2);

// Predict the value of y when x == 5.5, using a robustness level of 2.
auto prediction = model.predict(5.5, 2);
---

References:

Cleveland, W.S. (1979). "Robust Locally Weighted Regression and Smoothing
Scatterplots". Journal of the American Statistical Association 74 (368):
829-836.
*/
Loess1D loess1D(RX, RY)(
    RY y,
    RX x,
    double span,
    int degree = 1
) {
    dstatsEnforce(span > 0 && span <= 1, format(
        "Span must be >0 and <= 1 for loess1D, not %s.", span
    ));
    dstatsEnforce(degree >= 0, "Degree must be >= 0 for loess1D.");

    auto ret = new Loess1D;
    ret._x = array(map!(to!double)(x));
    ret._y = array(map!(to!double)(y));
    qsort(ret._x, ret._y);

    ret.nNeighbors = to!int(ret.x.length * span);
    dstatsEnforce(ret.nNeighbors > degree, format(
        "Cannot do degree %s loess with a window of only %s points.  " ~
        "Increase the span parameter.", degree, ret.nNeighbors
    ));

    ret._degree = degree;
    auto xPoly = new double[][degree];
    if(degree > 0) xPoly[0] = ret._x;

    foreach(d; 2..degree + 1) {
        xPoly[d - 1] = ret._x.dup;
        xPoly[d - 1][] ^^= d;
    }

    ret.xPoly = xPoly;
    return ret;
}

unittest {
    auto x = [1, 2, 3, 4, 5, 6, 7, 8];
    auto y = [3, 2, 8, 2, 6, 9, 0, 1];

    auto loess1 = loess1D(y, x, 0.75, 1);

    // Values from R's lowess() function.  This gets slightly different
    // results than loess(), probably due to disagreements bout windowing
    // details.
    assert(approxEqual(loess1.predictions(0),
        [2.9193046, 3.6620295, 4.2229953, 5.2642335, 5.3433985, 4.4225636,
         2.7719778, 0.6643268]
    ));

    loess1 = loess1D(y, x, 0.5, 1);
    assert(approxEqual(loess1.predictions(0),
        [2.1615941, 4.0041736, 4.5642738, 4.8631052, 5.7136895, 5.5642738,
         2.8631052, -0.1977227]
    ));

    assert(approxEqual(loess1.predictions(2),
        [2.2079526, 3.9809030, 4.4752888, 4.8849727, 5.7260333, 5.4465225,
         2.8769120, -0.1116018]
    ));

    // Test 0th and 2nd order using R's loess() function since lowess() doesn't
    // support anything besides first degree.
    auto loess0 = loess1D(y, x, 0.5, 0);
    assert(approxEqual(loess0.predictions(0),
        [3.378961, 4.004174, 4.564274, 4.863105, 5.713689, 5.564274, 2.863105,
         1.845369]
    ));

    // Not testing the last point.  R's loess() consistently gets slightly
    // different answers for the last point than either this function or
    // R's lowess() for degree > 0.  (This function and R's lowess() agree
    // when this happens.)  It's not clear which is right but the differences
    // are small and not practically important.
    auto loess2 = loess1D(y, x, 0.75, 2);
    assert(approxEqual(loess2.predictions(0)[0..$ - 1],
        [2.4029984, 4.1021339, 4.8288941, 4.5523535, 6.0000000, 6.4476465,
         3.7669741]
    ));
}

/**
This class is returned from the loess1D function and holds the state of a
loess regression with one predictor variable.
*/
final class Loess1D {
private:
    double[] _x;
    double[][] xPoly;
    double[] _y;
    double[][] yHat;
    double[][] residualWeights;

    int _degree;
    int nNeighbors;

    size_t nearestNeighborX(double point) const {
        auto sortedX = assumeSorted(x);
        auto trisected = sortedX.trisect(point);

        if(trisected[1].length) return trisected[0].length;
        if(!trisected[0].length) return 0;
        if(!trisected[2].length) return trisected[0].length - 1;

        if(point - trisected[0][trisected[0].length - 1] <
           trisected[2][0] - point
        ) {
            return trisected[0].length - 1;
        }

        return trisected[0].length;
    }

    static void computexTWx(
        const double[][] xPoly,
        const(double)[] weights,
        ref DoubleMatrix covMatrix
    ) {
        foreach(i; 0..xPoly.length) foreach(j; 0..i + 1) {
            covMatrix[i + 1, j + 1] = covMatrix[j + 1, i + 1] = threeDot(
                xPoly[i], weights, xPoly[j]
            );
        }

        // Handle intercept terms
        foreach(i; 1..covMatrix.rows) {
            covMatrix[0, i] = covMatrix[i, 0] = dotProduct(weights, xPoly[i - 1]);
        }

        covMatrix[0, 0] = sum(weights);
    }

    static void computexTWy(
        const(double)[][] xPoly,
        const(double)[] y,
        const(double)[] weights,
        double[] ans
    ) {
        foreach(i; 0..xPoly.length) {
            ans[i + 1] = threeDot(xPoly[i], weights, y);
        }

        // Intercept:
        ans[0] = dotProduct(weights, y);
    }

public:
    /**
    Predict the value of y when x == point, using robustness iterations
    of the biweight procedure outlined in the reference to make the
    estimates more robust.

    Notes:

    This function is computationally intensive but may be called
    from multiple threads simultaneously.  When predicting a
    large number of points, a parallel foreach loop may be used.

    Before calling this function with robustness > 0,
    computeRobustWeights() must be called.  See this function for details.

    Returns:  The predicted y value.
    */
    double predict(double point, int robustness = 0) const {
        dstatsEnforce(cast(int) residualWeights.length >= robustness,
            "For robust loess1D estimates, computeRobustWeights() " ~
            "must be called with the proper robustness level before " ~
            "calling predict()."
        );

        auto alloc = newRegionAllocator();
        auto covMatrix = doubleMatrix(degree + 1, degree + 1, alloc);
        auto xTy = alloc.uninitializedArray!(double[])(degree + 1);
        auto betas = alloc.uninitializedArray!(double[])(degree + 1);
        auto xPolyNeighbors = alloc.array(xPoly);

        size_t firstNeighbor = nearestNeighborX(point);
        size_t lastNeighbor = firstNeighbor;

        while(lastNeighbor - firstNeighbor + 1 < nNeighbors) {
            immutable upperDiff = (lastNeighbor + 1 >= x.length) ?
                double.infinity : (x[lastNeighbor + 1] - point);
            immutable lowerDiff = (firstNeighbor == 0) ?
                double.infinity : (point - x[firstNeighbor - 1]);

            if(upperDiff < lowerDiff) {
                lastNeighbor++;
            } else {
                firstNeighbor--;
            }
        }

        foreach(j, col; xPoly) {
            xPolyNeighbors[j] = xPoly[j][firstNeighbor..lastNeighbor + 1];
        }
        auto yNeighbors = y[firstNeighbor..lastNeighbor + 1];
        immutable maxDist = computeMaxDist(x[firstNeighbor..lastNeighbor + 1], point);
        auto w = alloc.uninitializedArray!(double[])(yNeighbors.length);
        foreach(j, neighbor; x[firstNeighbor..lastNeighbor + 1]) {
            immutable diff = abs(point - neighbor);
            w[j] = max(0, (1 - (diff / maxDist) ^^ 3) ^^ 3);
            if(robustness > 0) {
                w[j] *= residualWeights[robustness - 1][j + firstNeighbor];
            }
        }

        computexTWx(xPolyNeighbors, w, covMatrix);
        computexTWy(xPolyNeighbors, w, yNeighbors, xTy);
        choleskySolve(covMatrix, xTy, betas);

        double ret = 0;
        foreach(d; 0..degree + 1) {
            ret += betas[d] * point ^^ d;
        }

        return ret;
    }

    /**
    Compute the weights for robust loess, for all robustness levels <= the
    robustness parameter.  This computation is embarrassingly parallel, so if a
    TaskPool is provided it will be parallelized.

    This function must be called before calling predict() with a robustness
    value > 0.  computeRobustWeights() must be called with a robustness level
    >= the robustness level predict() is to be called with.  This is not
    handled implicitly because computeRobustWeights() is very computationally
    intensive and because it modifies the state of the Loess1D object, while
    predict() is const.  Forcing computeRobustWeights() to be called explicitly
    allows multiple instances of predict() to be evaluated in parallel.
    */
    void computeRobustWeights(int robustness, TaskPool pool = null) {
        dstatsEnforce(robustness >= 0, "Robustness cannot be <0.");

        if(cast(int) yHat.length < robustness) {
            computeRobustWeights(robustness - 1, pool);
        }

        if(yHat.length >= robustness + 1) {
            return;
        }

        if(robustness > 0) {
            auto resid = new double[y.length];
            residualWeights ~= resid;
            resid[] = y[] - yHat[$ - 1][];

            foreach(ref r; resid) r = abs(r);
            immutable sixMed = 6 * median(resid);

            foreach(ref r; resid) {
                r = biweight(r / sixMed);
            }
        }

        yHat ~= new double[y.length];

        if(pool is null) {
            foreach(i, point; x) {
                yHat[robustness][i] = predict(point, robustness);
            }
        } else {
            foreach(i, point; pool.parallel(x, 1)) {
                yHat[robustness][i] = predict(point, robustness);
            }
        }
    }

    /**
    Obtain smoothed predictions of y at the values of x provided on creation
    of this object, for the given level of robustness.  Evaluating
    these is computationally expensive and may be parallelized by
    providing a TaskPool object.
    */
    const(double)[] predictions(int robustness, TaskPool pool = null) {
        dstatsEnforce(robustness >= 0, "Robustness cannot be <0.");
        computeRobustWeights(robustness, pool);
        return yHat[robustness];
    }

const pure nothrow @property @safe:

    /// The polynomial degree for the local regression.
    int degree() {
        return _degree;
    }

    /// The x values provided on object creation.
    const(double)[] x() {
        return _x;
    }

    /// The y values provided on object creation.
    const(double)[] y() {
        return _y;
    }
}

private:
double biweight(double x) pure nothrow @safe {
    if(abs(x) >= 1) return 0;
    return (1 - x ^^ 2) ^^ 2;
}

double computeMaxDist(const double[] stuff, double x) pure nothrow @safe {
    double ret = 0;
    foreach(elem; stuff) {
        ret = max(ret, abs(elem - x));
    }

    return ret;
}

double absMax(double a, double b) {
    return max(abs(a), abs(b));
}

LogisticRes logisticRegressImpl(T, V...)
(bool inference, double ridge, T yIn, V input) {
    auto alloc = newRegionAllocator();

    static if(isFloatingPoint!(V[$ - 1])) {
        alias input[$ - 1] conf;
        alias V[0..$ - 1] U;
        alias input[0..$ - 1] xIn;
        enforceConfidence(conf);
    } else {
        alias V U;
        alias input xIn;
        enum conf = 0.95;
    }

    static assert(!isInfinite!T, "Can't do regression with infinite # of Y's.");
    static if(isRandomAccessRange!T) {
        alias yIn y;
    } else {
        auto y = toBools(yIn, alloc);
    }

    static if(U.length == 1 && isRoR!U) {
        static if(isForwardRange!U) {
            auto x = toRandomAccessRoR(y.length, xIn, alloc);
        } else {
            auto x = toRandomAccessRoR(y.length, alloc.array(xIn), alloc);
        }
    } else {
        auto xx = toRandomAccessTuple(xIn, alloc);
        auto x = xx.expand;
    }

    typeof(return) ret;
    ret.betas.length = x.length;
    if(inference) ret.stdErr.length = x.length;
    ret.logLikelihood = doMLENewton(ret.betas, ret.stdErr, ridge, y, x);

    if(!inference) return ret;

    static bool hasNaNs(R)(R range) {
        return !filter!isNaN(range).empty;
    }

    if(isNaN(ret.logLikelihood) || hasNaNs(ret.betas) || hasNaNs(ret.stdErr)) {
        // Then we didn't converge or our data was defective.
        return ret;
    }

    ret.nullLogLikelihood = .priorLikelihood(y);
    double lratio = ret.logLikelihood - ret.nullLogLikelihood;

    // Compensate for numerical fuzz.
    if(lratio < 0 && lratio > -1e-5) lratio = 0;
    if(lratio >= 0 && x.length >= 2) {
        ret.overallP = chiSquareCDFR(2 * lratio, x.length - 1);
    }

    ret.p.length = x.length;
    ret.lowerBound.length = x.length;
    ret.upperBound.length = x.length;
    immutable nDev = -invNormalCDF((1 - conf) / 2);
    foreach(i; 0..x.length) {
        ret.p[i] = 2 * normalCDF(-abs(ret.betas[i]) / ret.stdErr[i]);
        ret.lowerBound[i] = ret.betas[i] - nDev * ret.stdErr[i];
        ret.upperBound[i] = ret.betas[i] + nDev * ret.stdErr[i];
    }

    return ret;
}

private void logisticRegressPenalizedImpl(Y, X...)
(double[] betas, double lasso, double ridge, Y y, X xIn) {
    static if(isRoR!(X[0]) && X.length == 1) {
        alias xIn[0] x;
    } else {
        alias xIn x;
    }

    auto alloc = newRegionAllocator();
    auto ps = alloc.uninitializedArray!(double[])(y.length);
    betas[] = 0;
    auto betasRaw = alloc.array(betas[1..$]);

    double oldLikelihood = -double.infinity;
    double oldPenalty2 = double.infinity;
    double oldPenalty1 = double.infinity;
    enum eps = 1e-6;
    enum maxIter = 1000;

    auto weights = alloc.uninitializedArray!(double[])(y.length);
    auto z = alloc.uninitializedArray!(double[])(y.length);
    auto xMeans = alloc.uninitializedArray!(double[])(x.length);
    auto xSds = alloc.uninitializedArray!(double[])(x.length);
    double zMean = 0;
    auto xCenterScale = alloc.uninitializedArray!(double[][])(x.length);
    foreach(ref col; xCenterScale) col = alloc.uninitializedArray!(double[])(y.length);

    // Puts x in xCenterScale, with weighted mean subtracted and weighted
    // biased stdev divided.  Also standardizes z similarly.
    //
    // Returns:  true if successful, false if weightSum is so small that the
    //           algorithm has converged for all practical purposes.
    bool doCenterScale() {
        immutable weightSum = sum(weights);
        if(weightSum < eps) return false;

        xMeans[] = 0;
        zMean = 0;
        xSds[] = 0;

        // Do copying.
        foreach(i, col; x) {
            copy(take(col, y.length), xCenterScale[i]);
        }

        // Compute weighted means.
        foreach(j, col; xCenterScale) {
            foreach(i, w; weights) {
                xMeans[j] += w * col[i];
            }
        }

        foreach(i, w; weights) {
            zMean += w * z[i];
        }

        xMeans[] /= weightSum;
        zMean /= weightSum;
        z[] -= zMean;
        foreach(i, col; xCenterScale) col[] -= xMeans[i];

        // Compute biased stdevs.
        foreach(j, ref sd; xSds) {
            sd = sqrt(meanStdev(xCenterScale[j]).mse);
        }

        foreach(i, col; xCenterScale) col[] /= xSds[i];
        return true;
    }

    // Rescales the beta coefficients to undo the effects of standardizing.
    void rescaleBetas() {
        betas[0] = zMean;
        foreach(i, b; betasRaw) {
            betas[i + 1] = b / xSds[i];
            betas[0] -= betas[i + 1] * xMeans[i];
        }
    }

    foreach(iter; 0..maxIter) {
        evalPs(betas[0], ps, betas[1..$], x);
        immutable lh = logLikelihood(ps, y);
        immutable penalty2 = ridge * reduce!"a + b * b"(0.0, betas);
        immutable penalty1 = lasso * reduce!"a + (b < 0) ? -b : b"(0.0, betas);

        if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
        && (oldPenalty2 - penalty2) / (absMax(penalty2, oldPenalty2) + 0.1) < eps
        && (oldPenalty1 - penalty1) / (absMax(penalty1, oldPenalty1) + 0.1) < eps) {
            return;
        } else if(isNaN(lh) || isNaN(penalty2) || isNaN(penalty1)) {
            betas[] = double.nan;
            return;
        }

        oldPenalty2 = penalty2;
        oldPenalty1 = penalty1;
        oldLikelihood = lh;

        foreach(i, ref w; weights) {
            w = (ps[i] * (1 - ps[i]));
        }

        z[] = betas[0];
        foreach(i, col; x) {
            static if(is(typeof(col) : const(double)[])) {
                z[] += col[] * betas[i + 1];
            } else {
                foreach(j, ref elem; z) {
                    elem += col[j] * betas[i + 1];
                }
            }
        }

        foreach(i, w; weights) {
            if(w == 0){
                z[i] = 0;
            } else {
                immutable double yi = (y[i] == 0) ? 0.0 : 1.0;
                z[i] += (yi - ps[i]) / w;
            }
        }

        immutable centerScaleRes = doCenterScale();

        // If this is false then weightSum is so small that all probabilities
        // are for all practical purposes either 0 or 1.  We can declare
        // convergence and go home.
        if(!centerScaleRes) return;

        if(lasso > 0) {
            // Correct for different conventions in defining ridge params
            // so all functions get the same answer.
            immutable ridgeCorrected = ridge * 2.0;
            coordDescent(z, xCenterScale, betasRaw,
                lasso, ridgeCorrected, weights);
        } else {
            // Correct for different conventions in defining ridge params
            // so all functions get the same answer.
            immutable ridgeCorrected = ridge * 2.0;
            ridgeLargeP(z, xCenterScale, ridgeCorrected, betasRaw, weights);
        }

        rescaleBetas();
    }

    immutable lh = logLikelihood(ps, y);
    immutable penalty2 = ridge * reduce!"a + b * b"(0.0, betas);
    immutable penalty1 = lasso * reduce!"a + (b < 0) ? -b : b"(0.0, betas);

    if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
    && (oldPenalty2 - penalty2) / (absMax(penalty2, oldPenalty2) + 0.1) < eps
    && (oldPenalty1 - penalty1) / (absMax(penalty1, oldPenalty1) + 0.1) < eps) {
        return;
    } else {
        // If we got here, we haven't converged.  Return NaNs instead of bogus
        // values.
        betas[] = double.nan;
    }
}

// Calculate the mean squared error of all ranges.  This is delicate, though,
// because some may be infinite and we want to stop at the shortest range.
double[] calculateMSEs(U...)(U xIn, RegionAllocator alloc) {
    static if(isRoR!(U[0]) && U.length == 1) {
        alias xIn[0] x;
    } else {
        alias xIn x;
    }

    size_t minLen = size_t.max;
    foreach(r; x) {
        static if(!isInfinite!(typeof(r))) {
            static assert(hasLength!(typeof(r)),
                "Ranges passed to doMLENewton should be random access, meaning " ~
                "either infinite or with length.");

            minLen = min(minLen, r.length);
        }
    }

    dstatsEnforce(minLen < size_t.max,
        "Can't do logistic regression if all of the ranges are infinite.");

    auto ret = alloc.uninitializedArray!(double[])(x.length);
    foreach(ti, range; x) {
        //Workaround for bug https://issues.dlang.org/show_bug.cgi?id=13151
        //can be removed once no more need for 2.066.* support
        static if(is(typeof(range.save) R == Take!T, T)) {
            auto saved = range.save;
            ret[ti] = meanStdev(R(saved.source, min(minLen, saved.maxLength))).mse;
        } else {
            ret[ti] = meanStdev(take(range.save, minLen)).mse;
        }
    }

    return ret;
}

double doMLENewton(T, U...)
(double[] beta, double[] stdError, double ridge, T y, U xIn) {
    // This big, disgusting function uses the Newton-Raphson method as outlined
    // in http://socserv.mcmaster.ca/jfox/Courses/UCLA/logistic-regression-notes.pdf
    //
    // The matrix operations are kind of obfuscated because they're written
    // using very low-level primitives and with as little temp space as
    // possible used.
    static if(isRoR!(U[0]) && U.length == 1) {
        alias xIn[0] x;
    } else {
        alias xIn x;
    }

    auto alloc = newRegionAllocator();

    double[] mses;  // Used for ridge.
    if(ridge > 0) {
        mses = calculateMSEs(x, alloc);
    }

    beta[] = 0;
    if(stdError.length) stdError[] = double.nan;

    auto ps = alloc.uninitializedArray!(double[])(y.length);

    double getPenalty() {
        if(ridge == 0) return 0;

        double ret = 0;
        foreach(i, b; beta) {
            ret += ridge * mses[i] * b ^^ 2;
        }

        return ret;
    }

    enum eps = 1e-6;
    enum maxIter = 1000;

    auto oldLikelihood = -double.infinity;
    double oldPenalty = -double.infinity;
    auto firstDerivTerms = alloc.uninitializedArray!(double[])(beta.length);

    // matSaved saves mat for inverting to find std. errors, only if we
    // care about std. errors.
    auto mat = doubleMatrix(beta.length, beta.length, alloc);
    DoubleMatrix matSaved;

    if(stdError.length) {
        matSaved = doubleMatrix(beta.length, beta.length, alloc);
    }

    void saveMat() {
        foreach(i; 0..mat.rows) foreach(j; 0..mat.columns) {
            matSaved[i, j] = mat[i, j];
        }
    }

    void doStdErrs() {
        if(stdError.length) {
            // Here, we actually need to invert the information matrix.
            // We can use mat as scratch space since we don't need it
            // anymore.
            invert(matSaved, mat);

            foreach(i; 0..beta.length) {
                stdError[i] = sqrt(mat[i, i]);
            }
        }
    }

    auto updates = alloc.uninitializedArray!(double[])(beta.length);
    foreach(iter; 0..maxIter) {
        evalPs(ps, beta, x);
        immutable lh = logLikelihood(ps, y);
        immutable penalty = getPenalty();

        if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
        && (oldPenalty - penalty) / (absMax(penalty, oldPenalty) + 0.1) < eps) {
            doStdErrs();
            return lh;
        } else if(isNaN(lh)) {
            beta[] = double.nan;
            return lh;
        }

        oldLikelihood = lh;
        oldPenalty = penalty;

        foreach(i; 0..mat.rows) foreach(j; 0..mat.columns) {
            mat[i, j] = 0;
        }

        // Calculate X' * V * X in the notation of our reference.  Since
        // V is a diagonal matrix of ps[] * (1.0 - ps[]), we only have one
        // dimension representing it.  Do this for the lower half, then
        // symmetrize the matrix.
        foreach(i, xi; x) foreach(j, xj; x[0..i + 1]) {
            foreach(k; 0..ps.length) {
                mat[i, j] += (ps[k] * (1 - ps[k])) * xi[k] * xj[k];
            }
        }

        symmetrize(mat);

        if(stdError.length) {
            saveMat();
        }

        // Convert ps to ys - ps.
        foreach(pIndex, ref p; ps) {
            p = (y[pIndex] == 0) ? -p : (1 - p);
        }

        // Compute X'(y - p).
        foreach(ti, xRange; x) {
            //Workaround for bug https://issues.dlang.org/show_bug.cgi?id=13151
            //can be removed once no more need for 2.066.* support
            static if(is(typeof(xRange) R == Take!T, T)) {
                firstDerivTerms[ti] = dotProduct(
                        R(xRange.source, min(ps.length, xRange.maxLength)),
                        ps);
            } else {
                firstDerivTerms[ti] = dotProduct(take(xRange, ps.length), ps);
            }
        }

        // Add ridge penalties, if any.
        if(ridge > 0) {
            foreach(diagIndex, mse; mses) {
                mat[diagIndex, diagIndex] += 2 * ridge * mse;
                firstDerivTerms[diagIndex] -= beta[diagIndex] * 2 * ridge * mse;
            }
        }

        choleskySolve(mat, firstDerivTerms, updates);
        beta[] += updates[];

        debug(print) writeln("Iter:  ", iter);
    }

    immutable lh = logLikelihood(ps, y);
    immutable penalty = getPenalty();
    if((lh - oldLikelihood) / (absMax(lh, oldLikelihood) + 0.1) < eps
    && (oldPenalty - penalty) / (absMax(penalty, oldPenalty) + 0.1) < eps) {
        doStdErrs();
        return lh;
    } else {
        // If we got here, we haven't converged.  Return NaNs instead of bogus
        // values.
        beta[] = double.nan;
        return double.nan;
    }
}

private double logLikelihood(Y)(double[] ps, Y y) {
    double sum = 0;
    size_t i = 0;
    foreach(yVal; y) {
        scope(exit) i++;
        if(yVal) {
            sum += log(ps[i]);
        } else {
            sum += log(1 - ps[i]);
        }
    }
    return sum;
}

void evalPs(X...)(double[] ps, double[] beta, X xIn) {
    evalPs(0, ps, beta, xIn);
}

void evalPs(X...)(double interceptTerm, double[] ps, double[] beta, X xIn) {
    static if(isRoR!(X[0]) && X.length == 1) {
        alias xIn[0] x;
    } else {
        alias xIn x;
    }

    assert(x.length == beta.length);
    ps[] = interceptTerm;

    foreach(i, range; x) {
        static if(hasLength!(typeof(range))) {
            assert(range.length == ps.length);
        }

        static if(is(typeof(range) == double[])) {
            // Take advantage of array ops.
            ps[] += range[0..ps.length] * beta[i];
        } else {
            immutable b = beta[i];

            size_t j = 0;
            foreach(elem; range) {
                if(j >= ps.length) break;
                ps[j++] += b * elem;
            }
        }
    }

    foreach(ref elem; ps) elem = logistic(elem);
}

template isRoR(T) {
    static if(!isInputRange!T) {
        enum isRoR = false;
    } else {
        enum isRoR = isInputRange!(typeof(T.init.front()));
    }
}

template isFloatMat(T) {
    static if(is(T : const(float[][])) ||
        is(T : const(real[][])) || is(T : const(double[][]))) {
        enum isFloatMat = true;
    } else {
        enum isFloatMat = false;
    }
}

template NonRandomToArray(T) {
    static if(isRandomAccessRange!T) {
        alias T NonRandomToArray;
    } else {
        alias Unqual!(ElementType!(T))[] NonRandomToArray;
    }
}

double priorLikelihood(Y)(Y y) {
    uint nTrue, n;
    foreach(elem; y) {
        n++;
        if(elem) nTrue++;
    }

    immutable p = cast(double) nTrue / n;
    return nTrue * log(p) + (n - nTrue) * log(1 - p);
}

bool[] toBools(R)(R range, RegionAllocator alloc) {
    return alloc.array(map!"(a) ? true : false"(range));
}

auto toRandomAccessRoR(T)(size_t len, T ror, RegionAllocator alloc) {
    static assert(isRoR!T);
    alias ElementType!T E;
    static if(isArray!T && isRandomAccessRange!E) {
        return ror;
    } else static if(!isArray!T && isRandomAccessRange!E) {
        // Shallow copy so we know it has cheap slicing and stuff,
        // even if it is random access.
        return alloc.array(ror);
    } else {
        alias ElementType!E EE;
        auto ret = alloc.uninitializedArray!(EE[])(walkLength(ror.save));

        foreach(ref col; ret) {
            scope(exit) ror.popFront();
            col = alloc.uninitializedArray!(EE[])(len);

            size_t i;
            foreach(elem; ror.front) {
                col[i++] = elem;
            }
        }

        return ret;
    }
}

auto toRandomAccessTuple(T...)(T input, RegionAllocator alloc) {
    Tuple!(staticMap!(NonRandomToArray, T)) ret;

    foreach(ti, range; input) {
        static if(isRandomAccessRange!(typeof(range))) {
            ret.field[ti] = range;
        } else {
            ret.field[ti] = alloc.array(range);
        }
    }

    return ret;
}

// Borrowed and modified from Phobos's dotProduct() function,
// Copyright Andrei Alexandrescu 2008 - 2009.  Originally licensed
// under the Boost Software License 1.0.  (License text included
// at the beginning of this file.)
private double threeDot(
    const double[] x1,
    const double[] x2,
    const double[] x3
) in {
    assert(x1.length == x2.length);
    assert(x2.length == x3.length);
} body {
    immutable n = x1.length;
    auto avec = x1.ptr, bvec = x2.ptr, cvec = x3.ptr;
    typeof(return) sum0 = 0, sum1 = 0;

    const all_endp = avec + n;
    const smallblock_endp = avec + (n & ~3);
    const bigblock_endp = avec + (n & ~15);

    for (; avec != bigblock_endp; avec += 16, bvec += 16, cvec += 16)
    {
        sum0 += avec[0] * bvec[0] * cvec[0];
        sum1 += avec[1] * bvec[1] * cvec[1];
        sum0 += avec[2] * bvec[2] * cvec[2];
        sum1 += avec[3] * bvec[3] * cvec[3];
        sum0 += avec[4] * bvec[4] * cvec[4];
        sum1 += avec[5] * bvec[5] * cvec[5];
        sum0 += avec[6] * bvec[6] * cvec[6];
        sum1 += avec[7] * bvec[7] * cvec[7];
        sum0 += avec[8] * bvec[8] * cvec[8];
        sum1 += avec[9] * bvec[9] * cvec[9];
        sum0 += avec[10] * bvec[10] * cvec[10];
        sum1 += avec[11] * bvec[11] * cvec[11];
        sum0 += avec[12] * bvec[12] * cvec[12];
        sum1 += avec[13] * bvec[13] * cvec[13];
        sum0 += avec[14] * bvec[14] * cvec[14];
        sum1 += avec[15] * bvec[15] * cvec[15];
    }

    for (; avec != smallblock_endp; avec += 4, bvec += 4, cvec += 4) {
        sum0 += avec[0] * bvec[0] * cvec[0];
        sum1 += avec[1] * bvec[1] * cvec[1];
        sum0 += avec[2] * bvec[2] * cvec[2];
        sum1 += avec[3] * bvec[3] * cvec[3];
    }

    sum0 += sum1;

    /* Do trailing portion in naive loop. */
    while (avec != all_endp)
    {
        sum0 += *avec * *bvec * *cvec;
        ++avec;
        ++bvec;
        ++cvec;
    }

    return sum0;
}

unittest {
    auto a = [1.0,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21];
    auto b = a.dup;
    b[] *= 2;
    auto c = b.dup;
    c[] *= 3;
    immutable ans1 = threeDot(a, b, c);

    double ans2 = 0;
    foreach(i; 0..21) {
        ans2 += a[i] * b[i] * c[i];
    }

    assert(approxEqual(ans1, ans2));
}