/**This file contains all of the plots types in Plot2kill, which can be
 * drawn on screen using plot2kill.figure.Figure.
 *
 * Copyright (C) 2010-2012 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 plot2kill.plot;

import plot2kill.figure, plot2kill.util, plot2kill.hierarchical;
import std.random, std.typetuple;

version(dfl) {
    public import plot2kill.dflwrapper;
} else {
    public import plot2kill.gtkwrapper;
}

/**Abstract base class for all types of plot objects.*/
abstract class Plot {
protected:
    // Top of the plot.
    double upperLim = -double.infinity;

    // Bottom of the plot.
    double lowerLim = double.infinity;

    // Leftmost limit of the plot.
    double leftLim = double.infinity;

    // Rightmost limit of the plot.
    double rightLim = -double.infinity;

    string _legendText;
    
    // Hook in case the plot needs to reset its limits after seeing 
    // the figure it's about to be drawn to.
    void resetLimits(Figure figure) {}

package:
    PlotSize measureLegend(Font legendFont, FigureBase fig) {
        if(!legendText().length) return PlotSize(0, 0);
        auto ret = fig.measureText(legendText(), legendFont);
        ret.width += legendSymbolSize + legendSymbolTextSpace;
        ret.height = max(ret.height, legendSymbolSize);
        return PlotSize(ret.width, ret.height);
    }

public:
    // This should be package but for some reason package functions can't
    // be abstract.
    abstract void drawLegendSymbol(FigureBase fig, PlotRect where);

    /* Draw the plot on Figure using the rectangular area described by the
     * integer parameters.
     */
    abstract void drawPlot(Figure, double, double, double, double);

    /**
    Returns true if this plot type supports legends, otherwise false.
    */
    bool hasLegend() {
        return true;
    }

    /**Convenience method that instantiates a Figure object with this plot.
     * Useful for creating single-plot figures w/o a lot of boilerplate.
     *
     * Examples:
     * ---
     * auto hist = Histogram([1,2,3,4,5], 3).toFigure;
     * hist.showAsMain();
     * ---
     */
    Figure toFigure()  {
        return new Figure(this);
    }

    /**Instantiates a Figure object with this plot and also places the default
     * axes and tick labeling and title for the plot type, if any, on the
     * Figure.  If a plot type has no default labeling, simply forwards to
     * toFigure().
     */
    Figure toLabeledFigure()  {
        return toFigure;
    }

    /**The leftmost point on the plot.*/
    double leftMost()  {
        return leftLim;
    }

    /**The rightmost point on the plot.*/
    double rightMost()  {
        return rightLim;
    }

    /**The topmost point on the plot.*/
    double topMost()  {
        return upperLim;
    }

    /**The bottommost point on the plot.*/
    double bottomMost()  {
        return lowerLim;
    }

    ///
    string legendText()() {
        return _legendText;
    }

    ///
    This legendText(this This)(string newText) {
        _legendText = newText;
        return cast(This) this;
    }
}

/**A basic bar plot.*/
class BarPlot : Plot {
private:
    double[] _centers;
    double[] _heights;
    double[] _lowerErrors;
    double[] _upperErrors;
    BarPlot _stackOn;
    bool _outlineBar = true;

    // Return the bottom of the bar at a given index.  This is zero unless
    // we have a stacked plot.
    double bottom(size_t index) pure nothrow {
        return (_stackOn is null) ? 0 :
            _stackOn.bottom(index) + _stackOn._heights[index];
    }

    void fixBounds() {
        // We should never have both error bars and stacking.  The API doesn't
        // allow this, but double-check here.
        enforce(_stackOn is null ||
            (_lowerErrors.length == 0 && upperErrors.length == 0));

        this.leftLim = reduce!min(double.infinity, this.centers) - width / 2;
        this.rightLim = reduce!max(-double.infinity, this.centers) + width / 2;
        lowerLim = double.infinity;
        upperLim = -double.infinity;

        foreach(i, height; _heights) {
            if(lowerErrors.length > 0) {
                lowerLim = min(height - lowerErrors[i], lowerLim);
            } else {
                lowerLim = min(height + bottom(i), lowerLim);
            }

            if(upperErrors.length > 0) {
                upperLim = max(height + upperErrors[i], upperLim);
            } else {
                upperLim = max(height + bottom(i), upperLim);
            }
        }

        // Don't blend error bars into axes.  Handle the case where a boundary
        // is 0 as a special case, since it would look worse to have the axis
        // not be exactly zero if all bars are positive or all are negative than
        // to blend an error bar into an axis.
        if(upperErrors.length || lowerErrors.length) {
            immutable pad = 0.01 * (upperLim - lowerLim);
            if(upperErrors.length && upperLim != 0) {
                upperLim += pad;
            }
            if(lowerErrors.length && lowerLim != 0) {
                lowerLim -= pad;
            }
        }

        if(lowerLim > 0) {
            lowerLim = 0;
        }

        if(upperLim < 0) {
            upperLim = 0;
        }
    }


    this(double[] centers, double[] heights, double width) {
        this._centers = centers;
        this._heights = heights;
        this.width = width;
        _barColor = getColor(0, 0, 255);
        fixBounds();
    }

protected:
    override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {

        mixin(toPixels);
        mixin(drawErrorMixin);

        immutable multiplier = plotHeight / (this.upperLim - this.lowerLim);
        auto brush = form.getBrush(_barColor);
        scope(exit) doneWith(brush);

        auto blackPen = form.getPen(getColor(0, 0, 0));
        scope(exit) doneWith(blackPen);

        foreach(i, center; centers) {
            immutable zeroPoint = toPixelsY(bottom(i));
            immutable height = heights[i];
            immutable left = center - width / 2;
            immutable right = center + width / 2;
            immutable leftPixels = toPixelsX(left);
            immutable rightPixels = toPixelsX(right);
            immutable widthPixels = rightPixels - leftPixels;

            immutable heightPixels = roundTo!int(abs(height) * multiplier);

            immutable startAt = (height > 0) ?
                                zeroPoint - heightPixels :
                                zeroPoint;
            form.fillClippedRectangle(brush, leftPixels,
                startAt, widthPixels, heightPixels);

            if(_outlineBar) {
                form.drawClippedRectangle(blackPen, leftPixels,
                    startAt, widthPixels, heightPixels);
            }

            // Do error bars.
            if(lowerErrors.length) {
                drawErrorBar(blackPen,
                    center, height, height - lowerErrors[i], width / 2);
            }
            if(upperErrors.length) {
                drawErrorBar(blackPen,
                    center, height, height + upperErrors[i], width / 2);
            }
        }

        if(lowerLim < 0) {
            immutable zeroPoint = toPixelsY(0);
            auto pen = form.getPen(getColor(0, 0, 0), 2);
            // Draw line across figure at the zero point.
            form.drawClippedLine(pen,
                PlotPoint(toPixelsX(leftLim), zeroPoint),
                PlotPoint(toPixelsX(rightLim), zeroPoint)
            );
        }
    }

    override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        drawFillLegend(_barColor, fig, where);
    }

    Color _barColor;

public:
    /**Controls the color of the bar.  Defaults to blue.*/
    final Color barColor()() {
        return _barColor;
    }

    /// Setter
    final This barColor(this This)(Color newColor) {
        _barColor = newColor;
        return cast(This) this;
    }

    /**Controls whether each bar is outlined in a black rectange.*/
    final bool outlineBar()() {
        return _outlineBar;
    }

    final This outlineBar(this This)(bool shouldOutline) {
        _outlineBar = shouldOutline;
        return cast(This) this;
    }

    ///
    immutable double width;

    /**Create a BarPlot.  centers and heights must be input ranges with elements
     * implicitly convertible to double.  width determines the width of each
     * bar relative to the X-axis scale and must be greater than 0.
     */
    static BarPlot opCall(R1, R2)(R1 centers, R2 heights, double width)
    if(isInputRange!R1 && is(ElementType!R1 : double) &&
       isInputRange!R2 && is(ElementType!R2 : double)) {

        auto c = toDoubleArray(centers);
        auto h = toDoubleArray(heights);

        enforce(c.length == h.length,
            "Centers and heights must be same length for bar plot.");

        enforce(width > 0, "Width must be >0 for bar plot.");
        auto ret = new typeof(return)(c, h, width);
        return ret;
    }

    /**
    Create a BarPlot with error bars.  lowerErrors and upperErrors
     must be input ranges with elements implicitly convertible to double for
     error bars to be shown.  Any other value, such as null or 0, will result
     in no error bars being shown.  Therefore, to only show, for example,
     upper erros, simply pass in null or 0 for the lower errors.

     To draw symmetric error bars, simply pass in the same range for
     lowerErrors and upperErrors.  However, note that if you do this,
     the range will need to be a forward range, not an input range.
     */
     static BarPlot opCall(R1, R2, R3, R4)
     (R1 centers, R2 heights, double width, R3 lowerErrors, R4 upperErrors)
     if(isInputRange!R1 && is(ElementType!R1 : double) &&
        isInputRange!R2 && is(ElementType!R2 : double)) {

        auto ret = opCall(centers, heights, width);
        static if(isForwardRange!R3) {
            ret._lowerErrors = toDoubleArray(lowerErrors.save);
        } else static if(isInputRange!R3) {
            ret._lowerErrors = toDoubleArray(lowerErrors);
        }

        static if(isForwardRange!R4) {
            ret._upperErrors = toDoubleArray(upperErrors.save);
        } else static if(isInputRange!R4) {
            ret._upperErrors = toDoubleArray(upperErrors);
        }

        enforce(ret.upperErrors.length == 0 || ret.upperErrors.length ==
            ret.centers.length, "Length of upperErrors must equal number of bars.");
        enforce(ret.lowerErrors.length == 0 || ret.lowerErrors.length ==
            ret.centers.length, "Length of lowerErrors must equal number of bars.");

        ret.fixBounds();
        return ret;
    }

    /**
    Create a BarPlot that is to be stacked on top of this one.  For a
    convenience function for creating such a plot, see stackedBar().

    Warning:  Creating a cycle (two BarPlots mutually on top of each other)
              will cause infinite loops, stack overflows and otherwise bad
              things.
    */
    BarPlot stack(R)(R heights) {
        auto ret = new typeof(return)
            (_centers, toDoubleArray(heights), width);
        ret._stackOn = this;
        ret.fixBounds();
        return ret;
    }

    /**Scale this object by a factor of scaleFactor in the X direction.
     * This is useful for getting it onto the same scale as another plot.
     */
    void scaleCenters(double scaleFactor) {
        enforce(isFinite(scaleFactor), "Can't scale by infinity or NaN.");
        _centers[] *= scaleFactor;
        fixBounds();
    }

    /**Scale this object by a factor of scaleFactor in the Y direction.
     * This is useful for getting it onto the same scale as another plot.*/
    void scaleHeights(double scaleFactor) {
        enforce(isFinite(scaleFactor), "Can't scale by infinity or NaN.");
        _heights[] *= scaleFactor;
        fixBounds();
    }

    /**Shift this graph by shiftBy units in the X direction.
     * This is useful for getting it onto the same scale as another plot.
     */
    void shiftCenters(double shiftBy) {
        enforce(isFinite(shiftBy), "Can't shift by infinity or NaN.");
        _centers[] += shiftBy;
        fixBounds();
    }

    /**Get the centers of the bars.*/
    final const(double)[] centers()  {
        return _centers;
    }

    /**Get the heights of the bars.*/
    final const(double)[] heights()  {
        return _heights;
    }

    ///
    final const(double)[] lowerErrors()  {
        return _lowerErrors;
    }

    ///
    final const(double)[] upperErrors()  {
        return _upperErrors;
    }

    /**The default labeling includes each center receiving its own x tick
     * label if there are <= 10 bars on the graph.
     */
    override Figure toLabeledFigure()  {
        auto ret = toFigure;
        if(centers.length <= 10) {
            ret.xTickLabels(centers);
        }

        return ret;
    }
}

private template isRoR(T) {
    enum isRoR = isInputRange!T && isInputRange!(ElementType!T);
}

/**
Create an array of bar plots that, when inserted into a Figure, will effectively
become a grouped bar plot.

Parameters:

centers:
An input range of the overall center of each group of bars.

data:
A range of ranges, with one range for each bar color.  This range should have
one number for each group.

width:
The combined width of all of the bars for each group.

legendText:
An array of strings, one for each bar color, or null if no legend is desired.

colors:
An array of colors, one for each bar.  If none is provided, the default
colors are used.  These are, in order, blue, red, green, black, orange,
and purple.  If more colors are needed, they are generated using a random
number generator with a deterministic seed.

Examples:
---

// Make a plot with three groups:  One for "In Meeting", one for "On Phone",
// and one for "Coding".  The plot will also have two bar colors:  One for
// "Without Caffeine" and one for "With Caffeine".
auto withoutCaffeine = [8, 6, 3];
auto withCaffeine = [5, 3, 1];
auto sleepinessPlot = groupedBar(
    iota(3), [withoutCaffeine, withCaffeine], 0.6,
    ["W/o Caffeine", "W/ Caffeine"],
    [getColor(64, 64, 255), getColor(255, 64, 64)]
);
auto sleepinessFig = Figure(sleepinessPlot)
    .title("Sleepiness Survey")
    .yLabel("Sleepiness Rating")
    .xLabel("Activity")
    .legendLocation(LegendLocation.right)
    .horizontalGrid(true)
    .xTickLabels(
        iota(3),
        ["In Meeting", "On Phone", "Coding"]
    );
---
*/
BarPlot[] groupedBar(R1, R2)(
    R1 centers,
    R2 data,
    double width,
    string[] legendText = null,
    Color[] colors = null
)
if(isInputRange!R1 && isRoR!R2) {
     return groupedBar(centers, data, (double[][]).init, (double[][]).init,
        width, legendText, colors);
}

/**
Create a grouped bar plot with error bars.  lowerErrors and upperErrors
must either have the same dimensions as data or be empty.
*/
BarPlot[] groupedBar(R1, R2, R3, R4)(
    R1 centers,
    R2 data,
    R3 lowerErrors,
    R4 upperErrors,
    double width,
    string[] legendText = null,
    Color[] colors = null
)
if(isInputRange!R1 && isRoR!R2 && isRoR!R3 && isRoR!R4) {
    return groupStackImpl(CompoundBar.group, centers, data, lowerErrors,
        upperErrors, width, legendText, colors);
}

/**
Create a stacked bar plot.  The usage mechanics are identical to those of the
no error bar overload of groupedBar().

Examples:
---
// Stack coffee and tea consumption on top of each other.
Figure(
    stackedBar(iota(3), [[5, 3, 1], [1, 2, 3]], 0.6,
        ["Coffee", "Tea"]
    )
).legendLocation(LegendLocation.right)
    .title("Caffeine Consumption")
    .xLabel("Time of Day")
    .xTickLabels(iota(3), ["Morning", "Afternoon", "Evening"])
    .yLabel("Beverages")
    .showAsMain();
---
*/
BarPlot[] stackedBar(R1, R2)(
    R1 centers,
    R2 data,
    double width,
    string[] legendText = null,
    Color[] colors = null
)
if(isInputRange!R1 && isRoR!R2) {
    return groupStackImpl(CompoundBar.stack, centers, data, (double[][]).init,
        (double[][]).init, width, legendText, colors);
}

private enum CompoundBar {
    group,
    stack
}

private BarPlot[] groupStackImpl(R1, R2, R3, R4)(
    CompoundBar whatToDo,
    R1 centers,
    R2 data,
    R3 lowerErrors,
    R4 upperErrors,
    double width,
    string[] legendText = null,
    Color[] colors = null
) {
    if(whatToDo == CompoundBar.stack) {
        enforce(lowerErrors.length == 0 && upperErrors.length == 0);
    }

    auto centerArr = toDoubleArray(centers);
    auto dataArr = array(map!toDoubleArray(data));
    auto lerrArr = array(map!toDoubleArray(lowerErrors));
    auto uerrArr = array(map!toDoubleArray(upperErrors));

    foreach(elem; tuple(lerrArr, uerrArr, legendText, colors)) {
        enforce(elem.empty || elem.length == dataArr.length,
            "Range length mismatch in groupedBar.");
    }

    if(!colors.length) {
        colors = [
            getColor(0, 0, 255), getColor(255, 0, 0), getColor(0, 255, 0),
            getColor(0, 0, 0), getColor(255, 128, 0), getColor(255, 0, 255)
        ];

        // If we need more colors, generate them "randomly" but from a
        // deterministic seed.
        auto gen = Random(31415);
        while(colors.length < dataArr.length) {
            colors ~= getColor(
                uniform!"[]"(cast(ubyte) 0, cast(ubyte) 255, gen),
                uniform!"[]"(cast(ubyte) 0, cast(ubyte) 255, gen),
                uniform!"[]"(cast(ubyte) 0, cast(ubyte) 255, gen)
            );
        }
    }

    BarPlot[] ret;
    foreach(int groupIndex, group; dataArr) {
        BarPlot plot;

        enforce(group.length == centerArr.length,
            "Each group's length must be equal to centers.length for "
            ~ "grouped and stacked bar plots."
        );

        if(whatToDo == CompoundBar.group) {
            immutable groupWidth = width / dataArr.length;

            // Shift groupCenters over from overall centers.
            auto groupCenters = centerArr.dup;
            immutable offset = (groupIndex + 0.5) * groupWidth;
            groupCenters[] += offset - width / 2;

            if(lowerErrors.length && upperErrors.length) {
                plot = BarPlot(NoCopy(groupCenters), NoCopy(group), groupWidth,
                    NoCopy(lerrArr[groupIndex]), NoCopy(uerrArr[groupIndex])
                );
            } else {
                plot = BarPlot(NoCopy(groupCenters), NoCopy(group), groupWidth);
            }
        } else if(whatToDo == CompoundBar.stack) {
            if(groupIndex == 0) {
                plot = BarPlot(NoCopy(centerArr), NoCopy(group), width);
            } else {
                plot = ret[$ - 1].stack(NoCopy(group));
            }
        } else {
            assert(0);
        }

        plot.barColor(colors[groupIndex]);

        if(legendText.length) {
            plot.legendText(legendText[groupIndex]);
        }

        ret ~= plot;
    }

    return ret;
}

/**Determine behavior for elements outside of a fixed-border histogram's bounds.
 */
enum OutOfBounds {
     /** Throw throws an exception.*/
    throwException,

    /**Ignore simply skips the number.*/
    ignore,
    
    // For backwards compatibility, will eventually be removed.
    Throw = throwException,
    Ignore=ignore
}

/**Controls whether a histogram plots counts or probability.*/
enum HistType {
    /// The Y-axis should be counts.
    counts,

    /// The Y-axis should be probabilities.
    probability,
    
    // For backwards compatibility, will eventually be removed.
    Probability = probability,
    Counts = counts
}

/**A class for plotting regular (equal-width) histograms.*/
class Histogram : Plot {

    private double binWidth;

    private uint[] binCounts;
    private uint nElem;

    private HistType countsOrProbs;

    private this() {
        _barColor = getColor(0, 0, 255);
    }

    private bool isCumulative = false;

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        uint[] binCounts = this.binCounts;
        if(isCumulative) {
            binCounts = binCounts.dup;
            foreach(i; 1..binCounts.length) {
                binCounts[i] += binCounts[i - 1];
            }
        }

        immutable maxCount = reduce!max(0U, binCounts);
        immutable multiplier = plotHeight / cast(double) maxCount;

        immutable binWidth = cast(double) plotWidth / nBin;
        immutable bottom = plotHeight + topMargin;

        auto blackPen = form.getPen(getColor(0, 0, 0), 1);
        scope(exit) doneWith(blackPen);

        auto brush = form.getBrush(_barColor);
        scope(exit) doneWith(brush);

        double horizPos = leftMargin;
        double lastPosPixels = leftMargin;
        foreach(i, count; binCounts) {
            // Most of the complexity of this loop body is for making the bin
            // boundaries accurate in the context of having to round to
            // pixels.
            immutable barHeight = multiplier * count;
            immutable horizPixels = min(roundTo!int(horizPos), lastPosPixels);
            immutable stopAt = horizPos + binWidth;
            immutable thisBinWidth = max(1.0, stopAt - horizPixels);

            form.fillClippedRectangle(brush, lastPosPixels,
                bottom - barHeight, thisBinWidth, barHeight);
            form.drawClippedRectangle(blackPen, lastPosPixels,
                bottom - barHeight, thisBinWidth, barHeight);
            horizPos += binWidth;
            lastPosPixels = horizPixels + thisBinWidth;
        }
    }

    protected override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        drawFillLegend(_barColor, fig, where);
    }

    private void fixBounds() {
        if(isCumulative) {
            if(countsOrProbs == HistType.Probability) {
                upperLim = 1;
            } else {
                upperLim = nElem;
            }
        } else{
            if(countsOrProbs == HistType.Probability) {
                upperLim = reduce!max(0U, binCounts) / cast(double) nElem;
            } else {
                upperLim = reduce!max(0U, binCounts);
            }
        }
    }

    private OutOfBounds outOfBoundsBehavior = OutOfBounds.Throw;

    private Color _barColor;

    /**Controls the color of the bar.  Defaults to blue.*/
    final Color barColor()() {
        return _barColor;
    }

    /// Setter
    final This barColor(this This)(Color newColor) {
        _barColor = newColor;
        return cast(This) this;
    }


    /// The number of bins this histogram contains.
    final uint nBin() const pure nothrow  {
        return cast(uint) binCounts.length;
    }

    /**Factory method to instantiate this class.  nums must be a forward range
     * with elements implicitly convertible to double.  nBin specifies how many
     * bins the histogram should contain.
     */
    static Histogram opCall(R)(R nums, uint nBin)
    if(isForwardRange!R && is(ElementType!R : double)) {
        double leftLim = double.infinity, rightLim = -double.infinity;
        foreach(num; nums.save) {
            leftLim = min(leftLim, num);
            rightLim = max(rightLim, num);
        }

        return Histogram(nums, nBin, leftLim, rightLim);
    }

    /**Factory method to instantiate this class with predetermined limits.
     * This allows nums to be an input range instead of a forward range, since
     * no pass is necessary to compute the limits.
     *
     * This function both obeys and permanently sets whatever bounds behavior
     * is specified (throwing or ignoring on out of bounds numbers).  The
     * default behavior is to throw.  Rationale:  Errors should only pass
     * silently if explicitly silenced.
     */
    static Histogram opCall(R)(
        R nums,
        uint nBin,
        double leftLim,
        double rightLim,
        OutOfBounds outOfBoundsBehavior = OutOfBounds.throwException
    ) if(isInputRange!R && is(ElementType!R : double)) {
        auto ret = Histogram(nBin, leftLim, rightLim, outOfBoundsBehavior);

        foreach(num; nums) {
            ret.put(num);
        }

        return ret;
    }

    /**Create an empty histogram with pre-specified bounds, which will be
     * filled with data using the put method.
     *
     * Note:  The only reason this is a template is because of bugs in
     * overloading non-templated functions agsinst templated functions.
     */
    static Histogram opCall(I)(
        I nBin,
        double leftLim,
        double rightLim,
        OutOfBounds outOfBoundsBehavior = OutOfBounds.Throw
    ) if(isIntegral!I) {
        auto ret = new Histogram;
        ret.outOfBoundsBehavior = outOfBoundsBehavior;

        enforce(rightLim > leftLim,
            "Cannot create a histogram w/ upper lim <= lower lim.");
        enforce(nBin > 1, "Cannot create a histogram w/ <2 bins.");

        immutable binWidth = (rightLim - leftLim) / nBin;
        ret.binWidth = binWidth;

        ret.leftLim = leftLim;
        ret.rightLim = rightLim;
        ret.binCounts.length = to!uint(nBin);
        ret.lowerLim = 0;
        ret.upperLim = 0;

        return ret;
    }

    /**Add a number to the histogram on the fly.*/
    This put(this This)(double num) {
        if(outOfBoundsBehavior == OutOfBounds.Throw) {
            enforce(num >= leftLim && num <= rightLim, text(
                "Number out of bounds for histogram.  Got:  ", num,
                ", expected between ", leftLim, " and ", rightLim, "."));
        } else {
            if(!(num >= leftLim && num <= rightLim)) {
                return cast(This) this;
            }
        }

        uint bin;
        bin = to!uint((num - leftLim) / binWidth);
        if(bin == nBin) {  // Edge case.
            bin--;
        }

        binCounts[bin]++;
        nElem++;

        if(countsOrProbs == HistType.Counts) {
            if(isCumulative) {
                upperLim = nElem;
            } else if(binCounts[bin] > upperLim) {
                upperLim = binCounts[bin];
            }
        } else if(!isCumulative) {
            immutable binProb = binCounts[bin] / cast(double) nElem;
            if(binProb > upperLim) {
                upperLim = binProb;
            }
        }

        return cast(This) this;
    }

    /**Add the contents of another Histogram to this one.  The boundaries and
     * numbers of bins must be the same.  This histogram's settings are
     * retained.
     */
    This put(this This)(const Histogram rhs) {
        if(rhs is null) {
            return cast(This) this;
        }

        enforce(rhs.leftLim == this.leftLim && rhs.rightLim == this.rightLim,
            "Boundaries must be the same to combine histograms.");
        binCounts[] += rhs.binCounts[];
        nElem += rhs.nElem;

        fixBounds();
        return cast(This) this;
    }

    /**Assumes the LineGraph input is a plot of a PDF that this histogram is
     * supposed to approximate, and scales the Y axis of the LineGraph
     * accordingly so that both appear on the same scale.
     *
     * If this Histogram is cumulative, assumes that the input LineGraph is
     * a CDF instead.
     */
    This scaleDistributionFunction(this This)(LineGraph g) {
        if(isCumulative) {
            if(countsOrProbs == HistType.Counts) {
                g.scaleY(nElem);
            }
        // Don't need to do anything if this is a probability histogram.
        } else {
            double scaleFactor = (rightLim - leftLim) / nBin;
            if(countsOrProbs == HistType.Counts) {
                scaleFactor *= nElem;
            }

            g.scaleY(scaleFactor);
        }

        return cast(This) this;
    }

    /**Determine whether this object throws or ignores if it receives a number
     * outside its bounds via put.
     */
    This boundsBehavior(this This)(OutOfBounds behavior)  {
        outOfBoundsBehavior = behavior;
        return cast(This) this;
    }

    /**Set whether this histogram displays counts or probabilities.*/
    This histType(this This)(HistType newType) {
        countsOrProbs = newType;
        fixBounds();
        return cast(This) this;
    }

    /**Determines whether this histogram is cumulative.*/
    This cumulative(this This)(bool newVal) {
        isCumulative = newVal;
        fixBounds();
        return cast(This) this;
    }
}

/**Creates a histogram with equal frequency binning instead of equal width
 * binning.  The scale of a FrequencyHistogram is the probability density scale.
 *
 * Note that equal frequency binning doesn't work with discrete
 * distributions where probability density may be infinite at a point.
 * Therefore, if a probability density is calculated to be infinite, this
 * class will throw.
 *
 *
 */
class FrequencyHistogram : Plot {
    private double[] binWidths;
    private double elemsPerBin;

    private Color _barColor;

    this() {
        _barColor = getColor(0, 0, 255);
    }

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {

        mixin(toPixels);

        immutable multiplier = plotHeight / (this.upperLim - this.lowerLim);
        immutable zeroPoint = toPixelsY(0);
        auto brush = form.getBrush(_barColor);
        scope(exit) doneWith(brush);
        auto pen = form.getPen(_barColor, 1);
        scope(exit) doneWith(pen);

        auto blackPen = form.getPen(getColor(0, 0, 0));
        scope(exit) doneWith(blackPen);

        double xStart = leftLim;
        immutable totalWidth = rightLim - leftLim;
        foreach(width; binWidths) {
            scope(exit) xStart += width;
            immutable height = 1.0 / (binWidths.length * width);
            immutable leftPixels = toPixelsX(xStart);
            immutable rightPixels = toPixelsX(xStart + width);
            immutable widthPixels = rightPixels - leftPixels;
            immutable heightPixels = height * multiplier;

            immutable startAt = zeroPoint - heightPixels;

            form.fillClippedRectangle(brush, leftPixels,
                startAt, widthPixels, heightPixels);

            form.drawClippedRectangle(pen, leftPixels,
                startAt, widthPixels, heightPixels);

            // Don't outline rectangle because for equal-frequency
            // histograms this is more distracting than readable.
        }
    }

    override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        drawFillLegend(_barColor, fig, where);
    }

    /**Controls the color of the bar.  Defaults to blue.*/
    final Color barColor()() {
        return _barColor;
    }

    /// Setter
    final This barColor(this This)(Color newColor) {
        _barColor = newColor;
        return cast(This) this;
    }

    /**Create a FrequencyHistogram.  R must be an input range with elements
     * implicitly convertible to double.  nBin must be > 0 and <= nums.length.
     *
     * Throws:  Exception if the density for any bin is infinite.
     */
    static FrequencyHistogram opCall(R)(R nums, uint nBin) {
        auto numArr = toDoubleArray(nums);
        enforce(numArr.length >= nBin,
            "Can't create an equal-frequency histogram w/ < nBin elements.");

        sort(numArr);
        auto binWidths = new double[nBin];

        auto ret = new typeof(return);
        ret.leftLim = numArr[0];
        ret.rightLim = numArr[$ - 1];
        ret.lowerLim = 0;

        // Use linear interpolation to deal w/ non-integer elements/bin.
        immutable elemsPerBin = numArr.length / cast(double) nBin;
        double lastBinStop = numArr[0];
        foreach(i; 0..nBin) {
            immutable floatIndex = elemsPerBin * (i + 1);
            immutable indexFract = floatIndex - floor(floatIndex);
            size_t lowerIndex = to!size_t(floatIndex);
            size_t upperIndex = to!size_t(ceil(floatIndex));

            if(upperIndex >= numArr.length) {
                upperIndex = numArr.length - 1;
            }
            if(lowerIndex >= numArr.length) {
                lowerIndex = numArr.length - 1;
            }

            immutable diff = numArr[upperIndex] - numArr[lowerIndex];
            immutable curBinStop = numArr[lowerIndex] + diff * indexFract;

            binWidths[i] = curBinStop - lastBinStop;
            enforce(binWidths[i] > 0,
                "Can't create an equal-frequency histogram when some bin " ~
                "widths are 0."
            );
            lastBinStop = curBinStop;
        }

        ret.binWidths = binWidths;
        ret.elemsPerBin = elemsPerBin;
        immutable totalWidth = ret.rightLim - ret.leftLim;
        ret.upperLim = 1.0 / (binWidths.length * reduce!min(binWidths));
        return ret;
    }

    override Figure toLabeledFigure()  {
        auto ret = toFigure;
        ret.yLabel = "Probability Density";
        return ret;
    }
}

/**Creates a histogram in which every unique value gets its own bin.
 * Useful for histograms where the distribution is known to be discrete over
 * a small set of values.
 *
 * Hint:  Since this class inherits from BarPlot, BarPlot.centers will provide
 * a list of the unique values found.  This can be used to label the X
 * axis.
 */
class UniqueHistogram : BarPlot {

    /**The total count of this histogram.*/
    immutable uint nElem;

    private HistType countsOrProbs = HistType.counts;

    private this(double[] x, double[] y, double width, uint nElem) {
        this.nElem = nElem;
        super(x, y, width);
    }

    /**Create a UniqueHistogram.  R must be an input range with elements
     * implicitly convertible to double.  The width of each bin will be
     * widthFactor times the minimum distance between unique values. widthFactor
     * must be > 0 and <= 1.
     */
    static UniqueHistogram opCall(R)(R nums, double widthFactor = 0.8)
    if(isInputRange!R && isNumeric!(ElementType!R)) {
        enforce(widthFactor > 0 && widthFactor <= 1,
            "widthFactor must be > 0 and <= 1.");

        alias ElementType!R E;

        uint nElem;
        uint[double] counts;
        foreach(num; nums) {
            auto ptr = num in counts;
            if(ptr is null) {
                counts[num] = 1;
            } else {
                (*ptr)++;
            }
            nElem++;
        }

        auto possibleValues = counts.keys;
        sort(possibleValues);
        double minDiff = double.infinity;
        foreach(i; 1..possibleValues.length) {
            minDiff = min(minDiff, possibleValues[i] - possibleValues[i - 1]);
        }

        auto heights = new double[possibleValues.length];
        foreach(i, val; possibleValues) {
            heights[i] = counts[val];
        }

        return new typeof(this)
            (possibleValues, heights, minDiff * widthFactor, nElem);
    }

    /**Set whether this histogram displays counts or probabilities.*/
    void histType(HistType newType)  {
        if(newType == countsOrProbs) {
            return;
        } else if(newType == HistType.Counts) {
            scaleHeights(nElem);
        } else if(newType == HistType.Probability) {
            scaleHeights(1.0 / nElem);
        } else {
            assert(0);
        }

        countsOrProbs = newType;
    }  
}

/**Class for drawing a heat map.*/
class HeatMap : Plot {
    private double[][] values;
    private double minVal;
    private double maxVal;
    private uint _nRows;
    private uint _nCols;
    private Color[] _colors;

    enum string deprecatedMsg = "hotColor() and coldColor() are schedules " ~
        "for deprecation.  Use colors() instead.";

    private this() {
        // Set default colors.
        _colors = [
            getColor(255, 255, 255),
            getColor(0, 0, 0)
        ];
    }

    private Color getCellColor(double val)
    in {
        assert(_colors.length > 1);
        assert(val >= minVal);
        assert(val <= maxVal);
    } body {
        immutable diff = maxVal - minVal;
        immutable stride = diff / (_colors.length - 1.0);
        val -= minVal;
        val /= stride;

        immutable index1 = to!size_t(floor(val));
        immutable index2 = min(_colors.length - 1, to!size_t(ceil(val)));
        auto hotColor = _colors[index2];
        auto coldColor = _colors[index1];

        immutable fract = val - floor(val);
        immutable compl = 1.0 - fract;

        // Bug 4445:  roundTo!ubyte(255.0) throws.
        immutable red = cast(ubyte) roundTo!uint(
            coldColor.r * compl + hotColor.r * fract);
        immutable green = cast(ubyte) roundTo!uint(
            coldColor.g * compl + hotColor.g * fract);
        immutable blue = cast(ubyte) roundTo!uint(
            coldColor.b * compl + hotColor.b * fract);

        return getColor(red, green, blue);
    }

    protected void heatMapDefaultBounds() {
        enforceRectangular();

        _nRows = cast(uint) values.length;
        if(values.length > 0) {
            _nCols = cast(uint) values[0].length;
        }

        leftLim = 0.5;
        rightLim = nCols + 0.5;
        lowerLim = 0.5;
        upperLim = nRows + 0.5;
    }

    protected void setMinMax() {
        minVal = double.infinity;
        maxVal = -double.infinity;

        foreach(row; values) foreach(val; row) {
            minVal = min(minVal, val);
            maxVal = max(maxVal, val);
        }
    }

    protected override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        enforce(0, "Heat maps don't have legends.");
    }

    protected void enforceRectangular() {
        foreach(row; values) {
            enforce(row.length == values[0].length,
                "HeatMap matrices must be rectangular.");
        }
    }

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        immutable cellWidth = (rightLim - leftLim) / nCols;
        immutable cellHeight = (upperLim - lowerLim) / nRows;

        if(nRows < 1 || nCols < 1) {
            return;
        }

        mixin(toPixels);
        double lastRowStop = toPixelsY(upperLim);
        foreach(row; 0..nRows) {
            immutable rowStop = toPixelsY(upperLim - cellHeight * (row + 1));
            scope(exit) lastRowStop = rowStop;

            double lastColStop = toPixelsX(leftLim);
            foreach(col; 0..nCols) {
                immutable colStop = toPixelsX(cellWidth * (col + 1) + leftLim);
                scope(exit) lastColStop = colStop;

                auto color = getCellColor(values[row][col]);
                auto rect = PlotRect(lastColStop, lastRowStop,
                    colStop - lastColStop, rowStop - lastRowStop);

                auto brush = form.getBrush(color);
                scope(exit) doneWith(brush);
                auto pen = form.getPen(color, 1);
                scope(exit) doneWith(pen);

                form.fillClippedRectangle(brush, rect);
                form.drawClippedRectangle(pen, rect);
            }
        }
    }

    private Color _coldColor;
    private Color _hotColor;

    override bool hasLegend() {
        return false;
    }

    final Color coldColor()() {
        pragma(msg, deprecatedMsg);
        return _colors[0];
    }

    final This coldColor(this This)(Color newColor) {
        pragma(msg, deprecatedMsg);
        _colors[0] = newColor;
        return cast(This) this;
    }

    final Color hotColor()() {
        pragma(msg, deprecatedMsg);
        return _colors[1];
    }

    final This hotColor(this This)(Color newColor) {
        pragma(msg, deprecatedMsg);
        _colors[1] = newColor;
        return cast(This) this;
    }

    ///
    final Color[] colors()() {
        return _colors;
    }

    /**
    Set the colors to be used for this HeatMap as an array from coldest
    (smallest value) to warmest (largest value).  newColors.length must
    be >= 2.
    */
    final This colors(this This)(Color[] newColors) {
        enforce(newColors.length >= 2,
            "Need at least 2 colors for a HeatMap.");
        this._colors = newColors;
        return cast(This) this;
    }

    /**Create a heat map from a matrix represented as a range of ranges.
     * The matrix must be rectangular.  The elements of the ranges must
     * be implicitly convertible to double.
     */
    static HeatMap opCall(R)(R data)
    if(isInputRange!R && isInputRange!(ElementType!R) &&
       is(ElementType!(ElementType!(R)) : double)) {
        auto ret = new typeof(this);
        foreach(row; data) {
            ret.values ~= toDoubleArray(row);
        }

        ret.heatMapDefaultBounds();
        ret.setMinMax();
        return ret;
    }

    /**Create a heat map w/o copying data.*/
    static HeatMap noCopy(double[][] mat) {
        auto ret = new typeof(this);
        ret.values = mat;
        ret.heatMapDefaultBounds();
        ret.setMinMax();
        return ret;
    }

    ///
    final uint nRows()  {
        return _nRows;
    }

    ///
    final uint nCols()  {
        return _nCols;
    }

    /// Throws an exception.  Heat maps aren't allowed to have legends.
    override This legendText(this This)(string ignored) {
        enforce(0, "Heat maps can't have legend text.");
    }
}

/**
Convenience function that creates a hierarchically clustered heat map and,
if provided, rearranges your row and column names according to the clustering.

The distance and linkage aliases control the distance and linkage functions
for hierarchical clustering.  See plot2kill.hierarchical.hierarchicalCluster()
for details.
*/
HeatMap hierarchicalHeatMap(alias distance = euclidean, alias linkage = mean, R)
(R data, string[] rowNames = null, string[] colNames = null)
if(isInputRange!R && isInputRange!(ElementType!R) &&
is(ElementType!(ElementType!(R)) : double)) {

    auto matrix = array(map!(toDoubleArray)(data));
    enforce(matrix.length > 0 && matrix[0].length > 0,
        "Cannot produce a heat map w/ zero elements.");

    foreach(i; 0..matrix.length) {
        enforce(matrix[i].length == matrix[0].length,
            "data must be rectangular for hierarchicalHeatMap.");
    }

    auto rowClusters = hierarchicalCluster!(distance, linkage)
        (matrix, ClusterBy.rows, rowNames);

    auto rowPermApp = appender!(size_t[])();
    size_t nameIndex = 0;
    foreach(c; *rowClusters) {
        rowPermApp.put(c.index);
        if(rowNames.length) rowNames[nameIndex++] = c.name;
    }
    rowClusters = null;

    auto colClusters = hierarchicalCluster!(distance, linkage)
        (matrix, ClusterBy.columns, colNames);

    auto colPermApp = appender!(size_t[])();
    nameIndex = 0;
    foreach(c; *colClusters) {
        colPermApp.put(c.index);
        if(colNames.length) colNames[nameIndex++] = c.name;
    }
    colClusters = null;

    static T[] byPerm(T)(T[] input, const size_t[] perm) {
        assert(perm.length == input.length);

        auto ret = new T[input.length];
        foreach(i, p; perm) {
            ret[i] = input[p];
        }

        return ret;
    }

    assert(rowPermApp.data.length == matrix.length, text(
        rowPermApp.data.length, ' ', matrix.length));
    matrix = byPerm(matrix, rowPermApp.data);

    // This fixes some weirdness caused by different conventions between HeatMap
    // and most other plots.
    reverse(matrix);

    foreach(ref row; matrix) {
        row = byPerm(row, colPermApp.data);
    }

    return HeatMap.noCopy(matrix);
}

/**Creates a heat map representing the density of a 2-d probability
 * distribution.  This is useful when you want to visualize a joint
 * probability distribution but the sample size is so large that a
 * scatter plot would have an overwhelming number of points.
 */
class HeatScatter : HeatMap {
    private double cellWidth;
    private double cellHeight;

    private OutOfBounds outOfBoundsBehavior = OutOfBounds.Throw;

    /**Create a HeatScatter.  x, y must be forward ranges with elements
     * implicitly convertible to double, and must have the same lengths.
     */
    static
    HeatScatter opCall(R1, R2)(R1 x, R2 y, uint nRows = 10, uint nCols = 10)
    if(isForwardRange!R1 && isForwardRange!R2 &&
       is(ElementType!R1 : double) && is(ElementType!R2 : double)) {
        double xMin = double.infinity;
        double xMax = -double.infinity;
        uint xLen;
        foreach(num; x.save) {
            xLen++;
            xMin = min(num, xMin);
            xMax = max(num, xMax);
        }

        double yMin = double.infinity;
        double yMax = -double.infinity;
        uint yLen;
        foreach(num; y.save) {
            yLen++;
            yMin = min(num, yMin);
            yMax = max(num, yMax);
        }

        enforce(xLen == yLen,
            "Can't make HeatScatter when x.length != y.length.");
        return opCall(x, y, nRows, nCols, xMin, xMax, yMin, yMax);
    }

    /**Create a HeatScatter with pre-specified bounds.  x, y must be forward
     * ranges with elements implicitly convertible to double, and must have the
     * same lengths.
     */
    static
    HeatScatter opCall(R1, R2)(R1 x, R2 y, uint nRows, uint nCols,
    double xMin, double xMax, double yMin, double yMax,
    OutOfBounds boundsBehavior = OutOfBounds.Throw)
    if(isForwardRange!R1 && isForwardRange!R2 &&
       is(ElementType!R1 : double) && is(ElementType!R2 : double)) {

        auto ret = opCall(nRows, nCols, xMin, xMax, yMin, yMax, boundsBehavior);

        double maxVal = 0;
        while(!x.empty && !y.empty) {
            scope(exit) {
                x.popFront();
                y.popFront();
            }

            ret.put(x.front(), y.front());
        }

        enforce(x.empty && y.empty,
            "Can't make HeatScatter when x.length != y.length.");
        return ret;
    }

    /**Create a blank HeatScatter to fill in using the put() method.
     *
     * Note:  This is a template only because templates can't be overloaded
     * against regular functions.
     */
    static
    HeatScatter opCall(I)(I nRows, uint nCols,
        double xMin, double xMax, double yMin, double yMax,
        OutOfBounds boundsBehavior = OutOfBounds.Throw) {
        enforce(nRows > 0 && nCols > 0,
            "Can't make a heat scatter with 0 rows or columns.");

        enforce(xMax > xMin, "xMax must be > xMin.");
        enforce(yMax > yMin, "yMax must be > yMin.");

        auto grid = new double[][](nRows, nCols);
        foreach(row; grid) foreach(ref elem; row) {
            elem = 0;
        }

        auto ret = new typeof(return);
        ret.boundsBehavior = boundsBehavior;
        immutable cellWidth = (xMax - xMin) / nCols;
        immutable cellHeight = (yMax - yMin) / nRows;
        ret.cellWidth = cellWidth;
        ret.cellHeight = cellHeight;
        ret.values = grid;
        ret._nRows = nRows;
        ret._nCols = nCols;
        ret.minVal = 0;
        ret.maxVal = 0;
        ret.leftLim = xMin;
        ret.rightLim = xMax;
        ret.lowerLim = yMin;
        ret.upperLim = yMax;
        return ret;
    }

    /**Add an element to the plot.*/
    This put(this This)(double x, double y) {
        bool inBounds() {
            return (x >= leftLim && x <= rightLim && y >= lowerLim &&
            y <= upperLim);
        }

        if(outOfBoundsBehavior == OutOfBounds.Throw) {
            enforce(inBounds(), "Point out of bounds in HeatScatter.");
        } else if(!inBounds()) {
            return cast(This) this;
        }

        uint xCoord = to!uint((x - leftLim) / cellWidth);
        if(xCoord == nCols) {
            xCoord--;
        }

        uint yCoord = nRows - to!uint((y - lowerLim) / cellHeight) - 1;
        if(yCoord == uint.max) {
            yCoord = 0;
        }

        values[yCoord][xCoord]++;
        maxVal = max(maxVal, values[yCoord][xCoord]);

        return cast(This) this;
    }

    /**Add another HeatScatter's data to this.  The boundaries and row and
     * column counts must be the same.  The settings from this HeatScatter are
     * preserved.
     */
    This put(this This)(const HeatScatter rhs) {
        if(rhs is null) {
            return cast(This) this;
        }
        enforce(this.leftLim == rhs.leftLim && this.rightLim == rhs.rightLim &&
                this.upperLim == rhs.upperLim && this.lowerLim == rhs.lowerLim
                && this._nRows == rhs._nRows && this._nCols == rhs._nCols,
        "Cannot combine two HeatScatters w/ different bounds or row/column #s.");

        foreach(row; 0.._nRows) foreach(col; 0.._nCols) {
            this.values[row][col] += rhs.values[row][col];
            this.maxVal = max(this.values[row][col], this.maxVal);
        }

        return cast(This) this;
    }

    /**Determine whether this object throws or ignores if it receives a number
     * outside its bounds via put.
     */
    void boundsBehavior(OutOfBounds behavior) {
        outOfBoundsBehavior = behavior;
    }
}

/**Class for drawing a scatter plot.*/
class ScatterPlot : Plot {
    private double[] x;
    private double[] y;

    protected override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        drawTextLegend(_pointSymbol, _pointColor, fig, where);
    }

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        enforce(x.length == y.length);  // Should have already been checked.

        if(x.length < 1) {
            return;
        }

        mixin(toPixels);

        auto font = getFont(plot2kill.util.defaultFont,
            _pointSize + Figure.fontSizeAdjust);
        scope(exit) doneWith(font);

        auto pointDrawer = ScatterCharDrawer
            (_pointSymbol, font, _pointColor, form);

        pointDrawer.initialize;
        scope(exit) pointDrawer.restore();

        foreach(i; 0..x.length) {
            immutable curX = toPixelsX(x[i]);
            immutable curY = toPixelsY(y[i]);
            pointDrawer.draw(PlotPoint(curX, curY));
        }
    }

    private Color _pointColor;
    private char _pointSymbol = 'x';
    private int _pointSize = 10;

    /**The color of each point on the plot.*/
    final Color pointColor()() {
        return _pointColor;
    }

    /// Setter.
    final This pointColor(this This)(Color newColor) {
        _pointColor = newColor;
        return cast(This) this;
    }

    /**The symbol that should be used on the plot.  x and o work pretty well.
     * The default is x.
     */
    final char pointSymbol()() {
        return _pointSymbol;
    }

    /// Setter
    final This pointSymbol(this This)(char newSymbol) {
        _pointSymbol = newSymbol;
        return cast(This) this;
    }

    /// The size of a point.  (Default 10).
    final int pointSize()() {
        return _pointSize;
    }

    /// Setter
    final This pointSize(this This)(int newSize) {
        this._pointSize = newSize;
        return cast(This) this;
    }


    this() {
        _pointColor = getColor(0, 0, 0);
    }


    /**Factory method for creating a ScatterPlot.  x and y must both be
     * input ranges of the same length, with elements implicitly convertible
     * to doubles.  Note that they are copied inside the factory, so changes
     * to the original ranges after calling this factory will not affect the
     * plot.
     */
    static ScatterPlot opCall(R1, R2)(R1 x, R2 y)
    if(isInputRange!R1 && is(ElementType!R1 : double) &&
       isInputRange!R2 && is(ElementType!R2 : double)) {

        auto ret = new ScatterPlot;
        constructXYGraph(x, y, ret);

        fixXYGraphBounds(ret);
        addFudgeFactors(ret);
        return ret;
    }

    /**Convenience factory that produces a ScatterPlot with a default X
     * axis numbered 1, 2, ..., N where N is the number of points.  Mostly
     * useful for quick and dirty plots.
     */
    static ScatterPlot opCall(R)(R y)
    if(isInputRange!R && is(ElementType!R : double)) {
        auto ret = new ScatterPlot;
        ret.y = toDoubleArray(y);
        ret.x = new double[ret.y.length];

        foreach(i, ref elem; ret.x) {
            elem = i + 1;
        }

        fixXYGraphBounds(ret);
        addFudgeFactors(ret);
        return ret;
    }
}

/**Class for drawing a line graph, i.e. a set of points connected by lines.*/
class LineGraph : Plot {
    private double[] x;
    private double[] y;
    private double[] lowerErrors;
    private double[] upperErrors;


    private enum defaultErrorWidth = 0.05;
    private double _errorWidth = defaultErrorWidth;

    private char _pointSymbol = ' ';
    private Color _pointColor;
    private int _pointSize = 20;

    protected void fixBounds() {
        this.leftLim = reduce!min(double.infinity, this.x);
        this.rightLim = reduce!max(-double.infinity, this.x);

        this.upperLim = -double.infinity;
        this.lowerLim = double.infinity;
        foreach(i, yCoord; this.y) {
            if(this.lowerErrors.length) {
                this.lowerLim = min(this.lowerLim, yCoord - this.lowerErrors[i]);
            } else {
                this.lowerLim = min(this.lowerLim, yCoord);
            }
            if(this.upperErrors.length) {
                this.upperLim = max(this.upperLim, yCoord + this.upperErrors[i]);
            } else {
                this.upperLim = max(this.upperLim, yCoord);
            }
        }

        if(this.lowerErrors.length || this.upperErrors.length) {
            // Dont' cut off error bars.
            immutable yPad = (this.upperLim - this.lowerLim) * 0.01;
            immutable xPad = (this.rightLim - this.leftLim) * _errorWidth;
            this.leftLim -= xPad;
            this.rightLim += xPad;
            this.upperLim += yPad;
            this.lowerLim -= yPad;
        }

        if(_pointSymbol != ' ') {
 	        addFudgeFactors(this);
        }
    }

    protected override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        auto pen = fig.getPen(_lineColor, _lineWidth);
        scope(exit) doneWith(pen);
        drawLineLegend(pen, fig, where);

        if(_pointSymbol != ' ') {
            drawTextLegend(_pointSymbol, _pointColor, fig, where);
        }
    }

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        enforce(x.length == y.length);  // Should have already been checked.

        if(x.length < 2) {
            return;
        }

        mixin(toPixels);
        mixin(drawErrorMixin);

        immutable absErrorWidth = _errorWidth * (rightLim - leftLim);
        auto errorPen = form.getPen(getColor(0, 0, 0));
        scope(exit) doneWith(errorPen);

        double lastX = toPixelsX(x[0]);
        double lastY = toPixelsY(y[0]);

        void doErrors(size_t index) {
            if(lowerErrors.length) {
                drawErrorBar(
                    errorPen, x[index], y[index],
                        y[index] - lowerErrors[index], absErrorWidth
                );
            }
            if(upperErrors.length) {
                drawErrorBar(
                    errorPen, x[index], y[index],
                        y[index] + upperErrors[index], absErrorWidth
                );
            }
        }

        auto pen = form.getPen(_lineColor, _lineWidth);
        scope(exit) doneWith(pen);

        foreach(i; 1..x.length) {
            immutable curX = toPixelsX(x[i]);
            immutable curY = toPixelsY(y[i]);

            form.drawClippedLine(pen, PlotPoint(lastX, lastY), PlotPoint(curX, curY));
            lastX = curX;
            lastY = curY;
        }

        foreach(i; 0..x.length) {
            doErrors(i);
        }

        if(_pointSymbol == ' ') return;

        auto font = getFont(plot2kill.util.defaultFont,
            _pointSize + Figure.fontSizeAdjust);
        scope(exit) doneWith(font);

        auto pointDrawer = ScatterCharDrawer
            (_pointSymbol, font, _pointColor, form);

        pointDrawer.initialize;
        scope(exit) pointDrawer.restore();

        foreach(i; 0..x.length) {
            immutable curX = toPixelsX(x[i]);
            immutable curY = toPixelsY(y[i]);
            pointDrawer.draw(PlotPoint(curX, curY));
        }
    }

    private Color _lineColor;
    private uint _lineWidth = 1;

    /**The color of the line.  The default is black.*/
    final Color lineColor()() {
        return _lineColor;
    }

    /// Setter
    final This lineColor(this This)(Color newColor) {
        _lineColor = newColor;
        return cast(This) this;
    }

    /**The width of the line.  The default is 1.*/
    final uint lineWidth()() {
        return _lineWidth;
    }

    /// Setter
    final This lineWidth(this This)(uint newWidth) {
        _lineWidth = newWidth;
        return cast(This) this;
    }

    /**Error bar width, relative to the total width of the plot.  Must be
     * between 0 and 1.  If it's out of bounds, it will be set to the default
     * of 0.05.  If no error bars are to be drawn, this option is ignored.
     */
    final double errorWidth()() {
        return _errorWidth;
    }

    /// Setter
    final This errorWidth(this This)(double newWidth) {
        if(!(newWidth >= 0 && newWidth <= 1)) {
            _errorWidth = defaultErrorWidth;
        } else {
            _errorWidth = newWidth;
        }

        fixBounds();
        return cast(This) this;
    }

    /**
     The symbol that should be used to denote a data point.  Setting this
     to ' ' (a space) turns off point symbols and plots the line only.
     By default the point symbol is a space, so no point symbol is plotted..
     */
    final char pointSymbol()() {
        return _pointSymbol;
    }

    /// Setter
    final This pointSymbol(this This)(char newSymbol) {
        _pointSymbol = newSymbol;
        fixBounds();
        return cast(This) this;
    }

    /**The color of each point on the plot.*/
    final Color pointColor()() {
        return _pointColor;
    }

    /// Setter.
    final This pointColor(this This)(Color newColor) {
        _pointColor = newColor;
        return cast(This) this;
    }

    /// The size of a point.  (Default 20).
    final int pointSize()() {
        return _pointSize;
    }

    /// Setter
    final This pointSize(this This)(int newSize) {
        this._pointSize = newSize;
        fixBounds();
        return cast(This) this;
    }

    private this() {
        _lineColor = getColor(0, 0, 0);
        _pointColor = getColor(0, 0, 0);
    }

    /**Factory method for creating a LineGraph.  x and y must both be
     * input ranges of the same length, with elements implicitly convertible
     * to doubles.  Note that they are copied inside the factory, so changes
     * to the original ranges after calling this factory will not affect the
     * plot.
     */
    static LineGraph opCall(R1, R2)(R1 x, R2 y)
    if(isInputRange!R1 && is(ElementType!R1 : double) &&
       isInputRange!R2 && is(ElementType!R2 : double)) {

        auto ret = new LineGraph;
        constructXYGraph(x, y, ret);
        ret.fixBounds();
        return ret;
    }

    /**Convenience factory method that produces a LineGraph with a default
     * X axis numbered 1, 2, ..., N, where N is the number of data points,
     * and no error bars.  This is mostly useful for quick and dirty plots.
     */
    static LineGraph opCall(R)(R y)
    if(isInputRange!R && is(ElementType!R : double)) {
        auto ret = new LineGraph;
        ret.y = toDoubleArray(y);
        ret.x = new double[ret.y.length];

        foreach(i, ref elem; ret.x) {
            elem = i + 1;
        }
        ret.fixBounds();
        return ret;
    }

    /**Create a LineGraph with error bars.  lowerErrors and upperErrors
     * must be input ranges with elements implicitly convertible to double for
     * error bars to be shown.  Any other value, such as null or 0, will result
     * in no error bars being shown.  Therefore, to only show, for example,
     * upper erros, simply pass in null or 0 for the lower errors.
     *
     * To draw symmetric error bars, simply pass in the same range for
     * lowerErrors and upperErrors.  However, note that if you do this,
     * the range will need to be a forward range, not an input range.
     */
    static LineGraph opCall(R1, R2, R3, R4)
    (R1 x, R2 y, R3 lowerErrors, R4 upperErrors) {
        auto ret = new typeof(return);
        ret.x = toDoubleArray(x);
        ret.y = toDoubleArray(y);

        enforce(ret.x.length == ret.y.length,
            "x, y must have same length for line/scatter graph.");

        static if(isForwardRange!R3) {
            ret.lowerErrors = toDoubleArray(lowerErrors.save);
        } else static if(isInputRange!R3) {
            ret.lowerErrors = toDoubleArray(lowerErrors);
        }

        static if(isForwardRange!R4) {
            ret.upperErrors = toDoubleArray(upperErrors.save);
        } else static if(isInputRange!R4) {
            ret.upperErrors = toDoubleArray(upperErrors);
        }

        enforce(ret.upperErrors.length == 0 || ret.upperErrors.length ==
            ret.x.length, "Length of upperErrors must equal number of points.");
        enforce(ret.lowerErrors.length == 0 || ret.lowerErrors.length ==
            ret.x.length, "Length of lowerErrors must equal number of points.");

        if(ret.lowerErrors.length == 0 && ret.upperErrors.length == 0) {
            sort!"a[0] < b[0]"(zip(ret.x, ret.y));
        } else if(ret.lowerErrors.length == 0) {
            sort!"a[0] < b[0]"(zip(ret.x, ret.y, ret.upperErrors));
        } else if(ret.upperErrors.length == 0) {
            sort!"a[0] < b[0]"(zip(ret.x, ret.y, ret.lowerErrors));
        } else {
            sort!"a[0] < b[0]"(
                zip(ret.x, ret.y, ret.lowerErrors, ret.upperErrors));
        }

        ret.fixBounds();
        return ret;
    }

    /**Scale this object by a factor of scaleFactor in the X direction.
     * This is useful for getting it onto the same scale as another plot.
     */
    void scaleX(double scaleFactor) {
        enforce(isFinite(scaleFactor), "Can't scale by infinity or NaN.");
        x[] *= scaleFactor;
        fixBounds();
    }

    /**Scale this object by a factor of scaleFactor in the Y direction.
     * This is useful for getting it onto the same scale as another plot.*/
    void scaleY(double scaleFactor) {
        enforce(isFinite(scaleFactor), "Can't scale by infinity or NaN.");
        y[] *= scaleFactor;
        fixBounds();
    }

    /**Shift this graph by shiftBy units in the X direction.
     * This is useful for getting it onto the same scale as another plot.
     */
    void shiftX(double shiftBy) {
        enforce(isFinite(shiftBy), "Can't shift by infinity or NaN.");
        x[] += shiftBy;
        fixBounds();
    }

    /**Shift this graph by shiftBy units in the Y direction.
     * This is useful for getting it onto the same scale as another plot.
     */
    void shiftY(double shiftBy) {
        enforce(isFinite(shiftBy), "Can't shift by infinity or NaN.");
        y[] += shiftBy;
        fixBounds();
    }
}

/**
This class creates a best-fit line, and is useful for adding to scatter
plot figures.  The default line color is red instead of black and the default
width is two instead of one to make the best-fit line stand out against the
default scatter plot settings.
*/
class LinearFit : LineGraph {
    private double _alpha, _beta, _cor;
    
    final @property const pure nothrow @safe {
        /** 
        The intercept term of the fit line.
        */
        double alpha() { return _alpha; }
        
        /**
        The slope of the fit line.
        */
        double beta() { return _beta; }
        
        /**
        The correlation between the predictions made by the regression and
        the actual values.
        */
        double cor() { return _cor; }
        
        /**
        The R^2 value of the fit line.
        */
        double r2() { return _cor ^^ 2; }
    }

    private this() {
        super();
        _lineColor = getColor(255, 0, 0);
        _lineWidth = 2;
    }

    /**
    Compute the regression coefficients and create a best-fit line, where
    x is the independent variable and y is the dependent variable.
    lowerLim and upperLim control how far the best-fit line is extended in
    each direction along the X-axis.  By default the line is extended between 
    the minimum and maximum values of x.
    */
    static LinearFit opCall(R1, R2)
    (R1 x, R2 y, double lower = double.nan, double upper = double.nan)
    if(isInputRange!R1 && is(ElementType!R1 : double) &&
       isInputRange!R2 && is(ElementType!R2 : double)
    ) {
        
        static if(isForwardRange!R1) {
            alias x xx;
        } else {
            auto xx = toDoubleArray(x);
        }
        
        static if(isForwardRange!R2) {
            alias y yy;
        } else {
            auto yy = toDoubleArray(y);
        }
        
        import std.numeric;
        
        // This is somewhat inefficient, etc. because all the info has
        // to fit neatly on a scatter plot, so we can't be dealing with
        // millions of points, and because I don't want to add a dependency
        // to do it more properly.
        immutable xLen = walkLength(xx);
        immutable yLen = walkLength(yy);
        enforce(xLen == yLen, "x, y must be same length for LinearFit.");
        enforce(xLen > 1, "Need at least two elements for LinearFit.");
        
        immutable xMean = reduce!"a + b"(0.0, xx.save) / xLen;       
        immutable yMean = reduce!"a + b"(0.0, yy.save) / yLen;
    
        double xDotX = 0, xDotY = 0, yDotY = 0;
        foreach(xElem, yElem; lockstep(xx.save, yy.save)) {
            immutable xCenter = xElem - xMean;
            immutable yCenter = yElem - yMean;
            xDotX += xCenter ^^ 2;
            yDotY += yCenter ^^ 2;
            xDotY += xCenter * yCenter;
        }
        
        auto ret = new LinearFit();
        ret._beta = xDotY / xDotX;
        ret._alpha = yMean - ret._beta * xMean;
        
        ret._cor = xDotY / sqrt(xDotX) / sqrt(yDotY);
        
        if(isNaN(lower) || isNaN(upper)) {
            auto minMaxX = reduce!(min, max)(xx.save);
            lower = minMaxX[0];
            upper = minMaxX[1];
        }
        
        ret.x = [lower, upper];
        ret.y.length = 2;
        ret.y[] = ret.x[] * ret._beta + ret.alpha;
        ret.fixBounds();
        
        return ret;
    }
}

unittest {
    auto x = [8, 6, 7, 5, 3, 0, 9];    
    auto y = [3, 1, 4, 1, 5, 9, 2];
    
    // Make sure regression values are right.
    auto fit = LinearFit(x, y);
    assert(approxEqual(fit.cor, -0.7567996));
    assert(approxEqual(fit.beta, -0.6881));
    assert(approxEqual(fit.alpha, 7.3069));
    assert(approxEqual(fit.r2, 0.5727));
}

private void constructXYGraph(T, R1, R2)(R1 x, R2 y, T ret) {
    ret.x = toDoubleArray(x);
    ret.y = toDoubleArray(y);

    enforce(ret.x.length == ret.y.length,
        "x, y must have same length for line/scatter graph.");

    sort!"a[0] < b[0]"(zip(ret.x, ret.y));
    fixXYGraphBounds(ret);
}

private void fixXYGraphBounds(G)(G graph) {
    graph.lowerLim = reduce!min(double.infinity, graph.y);
    graph.upperLim = reduce!max(-double.infinity, graph.y);
    graph.leftLim = reduce!min(double.infinity, graph.x);
    graph.rightLim = reduce!max(-double.infinity, graph.x);

    /* Weird things happen when a graph has zero width or height.  Add fudge
     * factors to correct this.
     */
    with(graph) {
        if(lowerLim == upperLim) {
            lowerLim = nextafter(lowerLim, -double.infinity);
            upperLim = nextafter(upperLim, double.infinity);
        }

        if(leftLim == rightLim) {
            leftLim = nextafter(leftLim, -double.infinity);
            rightLim = nextafter(rightLim, double.infinity);
        }
    }
}

/**Plot a callable object on a range of values.*/
class ContinuousFunction : LineGraph {
    private enum string setupClass = q{
        auto ret = new ContinuousFunction;

        enforce(upperLim > lowerLim,
            "Can't do a function plot if upperLim <= lowerLim.");
        enforce(nEvals >= 2, "Can't do function plot with < 2 evals.");

        immutable diff = upperLim - lowerLim;
        immutable width = diff / (nEvals - 1);
        auto x = new double[nEvals];
        auto y = new double[nEvals];

        foreach(i; 0..nEvals) {
            x[i] = i * width + lowerLim;
            y[i] = callable(x[i]);
        }

        ret.x = x;
        ret.y = y;
        fixXYGraphBounds(ret);
        return ret;
    };

    /**Create a ContinuousFunction.  C is any callable type mapping a floating
     * point number to a number.  lowerLim is the lower limit of the plot.
     * upperLim is the upper limit of the plot.  nEvals is the number of
     * evalutations to perform.  The default is 1000.  More evaluations means
     * more accuracy but more computational intensity.
     */
    static ContinuousFunction opCall(C)(
        scope C callable,
        double lowerLim,
        double upperLim,
        uint nEvals = 1000
    ) if(is(typeof(C.init(2.0)) : double)) {
        // For some reason doing return opCall!callable doesn't work.
        // Using mixin instead.
        mixin(setupClass);
    }

    /**Create a ContinuousFunction using a template alias parameter instead of
     * a callable object.
     *
     * Note:  This function is given an name instead of opCall because DMD
     * thinks you're trying to instantiate a class if you do
     * ContinuousFunction!fun(...).
     */
    static ContinuousFunction fromAlias(alias callable)(
        double lowerLim,
        double upperLim,
        uint nEvals = 1000
    ) if(is(typeof(callable(2.0)) : double)) {
        mixin(setupClass);
    }

}

/**Plot a callable object on a range of values.*/
class DiscreteFunction : BarPlot {

    private this(double[] x, double[] y) {
        super(x, y, 1);
    }

    /**Create a DiscreteFunction. C is any callable type mapping an integer
     * to a number.  lowerLim is the lower limit of the plot.
     * upperLim is the upper limit of the plot.
     */
    static DiscreteFunction opCall(C)(
        scope C callable,
        int lowerLim,
        int upperLim,
    ) if(is(typeof(C.init(2)) : double)) {
        enforce(upperLim > lowerLim,
            "Can't do a function plot if upperLim <= lowerLim.");

        static if(is(typeof(C.opCall))) {
            alias ParameterTypeTuple!(C.opCall)[0] I;
        } else {
            alias ParameterTypeTuple!(C)[0] I;
        }

        I input = to!I(lowerLim);
        double[] x = new double[upperLim - lowerLim + 1];
        double[] y = new double[upperLim - lowerLim + 1];

        foreach(ulIndex; 0..x.length) {
            immutable index = cast(int) ulIndex;
            x[index] = lowerLim + index;
            y[index] = callable(input);
            input++;
        }

        auto ret = new typeof(return)(x, y);
        ret.fixBounds();
        return ret;
    }

    /**Create a DiscreteFunction using a template alias parameter instead of
     * a callable object.
     *
     * Note:  This function is given an name instead of opCall because DMD
     * thinks you're trying to instantiate a class if you do
     * ContinuousFunction!fun(...).
     */
    static DiscreteFunction fromAlias(alias callable)(
        int lowerLim,
        int upperLim,
    ) if(is(typeof(callable(2)) : double)) {
        enforce(upperLim > lowerLim,
            "Can't do a function plot if upperLim <= lowerLim.");

        int input = lowerLim;
        double[] x = new double[upperLim - lowerLim + 1];
        double[] y = new double[upperLim - lowerLim + 1];

        foreach(int index; 0..x.length) {
            x[index] = lowerLim + index;
            y[index] = callable(input);
            input++;
        }

        auto ret = new typeof(return)(x, y);
        ret.fixBounds();
        return ret;
    }
}

// Above was relatively basic plots.  Below is more specialized statistical
// plots.

/**Plots a ROC curve, or a curve with sensitivity on the Y-axis
 * and 1 - specificity on the X-axis.  This is a useful metric for
 * determining how well a test statistic discriminates between two classes.
 * The following assumptions are made in this implementation:
 *
 * 1.  For some cutoff value c and test statistic T, your decision rule is of
 *     the form "Class A if T larger than c, Class B if T smaller than c".
 *
 * 2.  In the case of ties, i.e. if class A and class B both have an identical
 *     value, linear interpolation is used.  This is because changing the
 *     value of c infinitesimally will change both sensitivity and specificity
 *     in these cases.
 */
class RocCurve : LineGraph {
private:
    double _auroc;

    this() {
        super();
    }

public:
    /**Create a RocCurve.  classATs are the test statistics that are
     * "supposed" to be bigger, classBTs are the test statistics that are
     * "supposed to be smaller.  Both R1 and R2 must be input ranges with
     * elements implicitly convertible to double.
     */
    static RocCurve opCall(R1, R2)(R1 classATs, R2 classBTs)
    if(isNumeric!(ElementType!R1) && isNumeric!(ElementType!R2)) {
        // Shamelessly cut-and-pasted and modified from dstats.

        auto classA = array(classATs);
        auto classB = array(classBTs);
        sort(classA);
        sort(classB);

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

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

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

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

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

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

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

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

            x ~= curPoint[0];
            y ~= curPoint[1];

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

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

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

        x ~= 1;
        y ~= 1;

        auto ret = new typeof(return)();
        ret._auroc = area;
        ret.x = x;
        ret.y = y;
        fixXYGraphBounds(ret);
        return ret;
    }

    /**Returns the area under the ROC curve.*/
    final double auroc() const pure nothrow  {
        return _auroc;
    }

    /**Default title is the area under the curve.  Default x label is
     * "1 - Specificity".  Default y label is "Sensitivity".
     */
    override Figure toLabeledFigure()  {
        auto ret = toFigure;
        ret.title = "AUROC = " ~ to!string(auroc);
        ret.xLabel = "1 - Specificity";
        ret.yLabel = "Sensitivity";
        return ret;
    }
}

unittest {
    // Values worked out by hand on paper.  If you don't believe me, work
    // them out yourself.
    auto foo = RocCurve([4,5,6], [1,2,3]);
    assert(foo.auroc == 1);

    foo = RocCurve([8,6,7,5,3,0,9], [3,6,2,4,3,6]);
    assert(approxEqual(foo.auroc, 0.6904762));

    foo = RocCurve([2,7,1,8,2,8,1,8], [3,1,4,1,5,9,2,6]);
    assert(approxEqual(foo.auroc, 0.546875));
}

/**Plots the quantiles of a set of data on the Y axis against the theoretical
 * qualtiles or the quantiles of another set of data on the X axis.
 */
class QQPlot : ScatterPlot {
    private this() {
        super();
        _lineColor = getColor(255, 0, 0);
    }

    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        super.drawPlot(
            form, leftMargin, topMargin, plotWidth, plotHeight
        );

        // Draw line that indicates identical distributions.
        if(max(lowerLim, leftLim) < min(upperLim, rightLim)) {
            mixin(toPixels);

            auto pen = form.getPen(_lineColor, _lineWidth);
            scope(exit) doneWith(pen);

            immutable lowerX = toPixelsX(max(lowerLim, leftLim));
            immutable lowerY = toPixelsY(max(lowerLim, leftLim));
            immutable upperX = toPixelsX(min(upperLim, rightLim));
            immutable upperY = toPixelsY(min(upperLim, rightLim));

            form.drawClippedLine(
                pen,
                PlotPoint(lowerX, lowerY),
                PlotPoint(upperX, upperY)
            );
        }
    }

    private Color _lineColor;

    /**The color of the y = x line that indicates identical distributions.
     * The default is red.
     */
    final Color lineColor()() {
        return _lineColor;
    }

    /// Setter
    final This lineColor(this This)(Color newColor) {
        _lineColor = newColor;
        return cast(This) this;
    }

    private double _lineWidth = 2;

    /**The width of the line that indicates identical distributions.
     * The default is 2.
     */
    final double lineWidth()() {
        return _lineWidth;
    }

    /// Setter.
    final This lineWidth(this This)(double newWidth) {
        _lineWidth = newWidth;
        return cast(This) this;
    }

    /**Create a QQPlot.  dataRange must be an input range with elements
     * implicitly convertible to doubles.  quantileFunction must be a
     * callable (function pointer, delegate, functor, etc.) mapping any
     * number between 0 and 1 to its quantile.
     *
     * Examples:
     * ---
     * auto norms = randArray!rNorm(100, 0, 1);
     * auto theoretical = paramFunctor!invNormalCDF(0, 1);
     * auto fig = new Figure(
     *     QQPlot(norms, theoretical)
     * );
     * fig.showAsMain();
     * ---
     */
    static QQPlot opCall(R, C)(R dataRange, C quantileFunction)
    if(is(typeof(C.init(2.0)) : double) && isInputRange!R &&
    is(ElementType!R : double)) {
        auto ret = new QQPlot;

        ret.y = toDoubleArray(dataRange);
        sort(ret.y);
        immutable double N = ret.y.length + 1.0;

        ret.x = new double[ret.y.length];
        foreach(i, ref elem; ret.x) {
            elem = quantileFunction((i + 1) / N);
        }

        fixXYGraphBounds(ret);
        addFudgeFactors(ret);
        return ret;
    }

    /**Default x label is "Theoretical Quantiles".  Default y label is
     * "Empirical Quantiles".
     */
    override Figure toLabeledFigure() {
        auto ret = toFigure;
        ret.xLabel = "Theoretical Quantiles";
        ret.yLabel = "Empirical Quantiles";
        return ret;
    }
}

/**
Draw a basic box-and-whisker plot.  The plots are drawn centered at Y
coordinates [0, 1, ..., N) unless the offset property is set to something else
(which is useful for putting multiple box pots with different color boxes on
the same figure).
*/
class BoxPlot : Plot {
private:
    double[] medians;
    double[] boxBottoms;
    double[] boxTops;
    double[] whiskerBottoms;
    double[] whiskerTops;
    double[][] outliers;
    double whiskerPercentile;
    double _offset = 0;
    Color _color;

    this() {
        _color = getColor(0, 0, 0);  // Default color.
    }

    void updateBounds() {
        if(medians.length == 0) return;

        upperLim = reduce!max(-double.infinity,
            chain(whiskerTops, join(outliers))
        );

        lowerLim = reduce!min(double.infinity,
            chain(whiskerBottoms, join(outliers))
        );

        immutable diff = upperLim - lowerLim;
        upperLim += 0.01 * diff;
        lowerLim -= 0.01 * diff;

        leftLim = -0.6;
        rightLim = medians.length - 0.4;

        leftLim += _offset;
        rightLim += _offset;

    }

    void addDataImpl(R)(R range) {
        auto doubles = toDoubleArray(range);
        enforce(doubles.length,
            "Cannot add a zero-length range to a box and whisker plot.");
        sort(doubles);

        if(doubles.length & 1) {
            medians ~= doubles[$ / 2];
        } else {
            medians ~= 0.5 * doubles[$ / 2] + 0.5 * doubles[$ / 2 + 1];
        }

        double doInterp(double percentile) {
            immutable floatIndex = percentile * (doubles.length - 1);
            immutable floored = to!size_t(floatIndex);
            immutable fract = floatIndex - floored;

            if(fract == 1 || floored == doubles.length - 1) {
                return doubles[floored];
            } else {
                return (1 - fract) * doubles[floored] +
                    fract * doubles[floored + 1];
            }
        }

        boxBottoms ~= doInterp(0.25);
        boxTops ~= doInterp(0.75);
        whiskerBottoms ~= doInterp(whiskerPercentile);
        whiskerTops ~= doInterp(1 - whiskerPercentile);

        immutable whiskerTopIndex = to!size_t(
            ceil((1 - whiskerPercentile) * (doubles.length - 1)) + 1
        );
        immutable whiskerBottomIndex = to!size_t(
            whiskerPercentile * (doubles.length - 1)
        );

        auto outlierRange = chain(
            doubles[0..whiskerBottomIndex],
            doubles[whiskerTopIndex..$]
        );

        outliers.length += 1;

        foreach(outlier; outlierRange) {
            outliers[$ - 1] ~= outlier;
        }

        // We're absolutely sure we don't need doubles anymore, and may
        // be working with huge datasets.
        delete doubles;

        updateBounds();
    }

protected:
    protected override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        mixin(toPixels);
        immutable nBoxes = medians.length;
        immutable boxWidth = 1;

        auto pen = form.getPen(_color, 1);
        auto font = getFont(plot2kill.util.defaultFont,
            10 + Figure.fontSizeAdjust);
        auto outlierDrawer = ScatterCharDrawer('o', font, _color, form);

        foreach(boxIndex, med; medians) {
            immutable left = boxWidth * boxIndex + 0.15 - 0.5 + _offset;
            immutable center = boxWidth * boxIndex + _offset;
            immutable right = boxWidth * boxIndex + 0.85 - 0.5 + _offset;
            immutable wLeft = boxWidth * boxIndex + 0.35 - 0.5 + _offset;
            immutable wRight = boxWidth * boxIndex + 0.65 - 0.5 + _offset;

            immutable leftPixels = toPixelsX(left);
            immutable centerPixels = toPixelsX(center);
            immutable rightPixels = toPixelsX(right);
            immutable wLeftPixels = toPixelsX(wLeft);
            immutable wRightPixels = toPixelsX(wRight);
            immutable boxTopPixels = toPixelsY(boxTops[boxIndex]);
            immutable boxBotPixels = toPixelsY(boxBottoms[boxIndex]);
            immutable wTopPixels = toPixelsY(whiskerTops[boxIndex]);
            immutable wBotPixels = toPixelsY(whiskerBottoms[boxIndex]);

            // Draw box.
            form.drawClippedRectangle(pen,
                leftPixels, boxTopPixels,
                rightPixels - leftPixels,
                boxBotPixels - boxTopPixels
            );

            // Draw median lines.
            form.drawClippedLine(
                pen,
                PlotPoint(leftPixels, toPixelsY(med)),
                PlotPoint(rightPixels, toPixelsY(med))
            );


            // Draw whiskers.
            form.drawClippedLine(
                pen,
                PlotPoint(centerPixels, boxTopPixels),
                PlotPoint(centerPixels, wTopPixels)
            );

            form.drawClippedLine(
                pen,
                PlotPoint(centerPixels, boxBotPixels),
                PlotPoint(centerPixels, wBotPixels)
            );

            form.drawClippedLine(
                pen,
                PlotPoint(wLeftPixels, wBotPixels),
                PlotPoint(wRightPixels, wBotPixels)
            );

            form.drawClippedLine(
                pen,
                PlotPoint(wLeftPixels, wTopPixels),
                PlotPoint(wRightPixels, wTopPixels)
            );

            if(whiskerPercentile == 0) continue;

            // Draw outliers.
            outlierDrawer.initialize();
            scope(exit) outlierDrawer.restore();

            foreach(outlier; outliers[boxIndex]) {
                outlierDrawer.draw(
                    PlotPoint(centerPixels, toPixelsY(outlier))
                );
            }
        }
    }

    override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        auto pen = fig.getPen(_color, 1);
        scope(exit) doneWith(pen);
        drawLineLegend(pen, fig, where);
    }

public:

    /**
    Create a BoxPlot.  whiskerPercentile controls the percentile at which
    data points are considered outliers and plotted as individual points.
    For any x, a percentile of x is equivalent to a percentile of 1 - x.
    For example, if whiskerPercentile is either 0 or 1, no data is plotted as
    individual points and the whiskers extend all the way to the extrema.
    */
    static BoxPlot opCall(double whiskerPercentile = 0) {
        enforce(whiskerPercentile >= 0 && whiskerPercentile <= 1,
            "whiskerPercentile must be between 0 and 1.");
        auto ret = new BoxPlot;
        ret.whiskerPercentile = min(whiskerPercentile, 1 - whiskerPercentile);
        return ret;
    }

    /**
    Add data to the box plot.  data may be any combination of ranges of
    ranges and individual ranges.
    */
    This addData(this This, R...)(R data)
    if(allSatisfy!(isInputRange, R)) {
        foreach(r; data) {
            static if(isInputRange!(ElementType!(typeof(r)))) {
                foreach(rr; r) addDataImpl(rr);
            } else {
                addDataImpl(r);
            }
        }

        return cast(This) this;
    }

    /// The offset from zero at which the first box is centered.
    double offset()() {
        return _offset;
    }

    /// Setter
    This offset(this This)(double newOffset) {
        this._offset = newOffset;
        updateBounds();
        return cast(This) this;
    }

    /// The color that the boxes, whiskers and outliers are drawn in.
    Color color()() {
        return _color;
    }

    /// Setter.
    This color(this This)(Color newColor) {
        this._color = newColor;
        return cast(This) this;
    }

    /// The number of boxes currently on the plot.
    int nBoxes() {
        // If you have more than a few billion boxes, you have bigger
        // problems than integer overflow.
        return cast(int) medians.length;
    }
}

/**
A dendrogram is a tree of hierarchically clustered data such that the height
of each split is proportional to the distance between clusters.
*/
class Dendrogram : Plot {
private:
    Cluster* clusters;
    immutable int nLeaves;
    string[] _names;
    
    // Make min leaf lines 10% of the difference between max, min distance.
    enum leafLineFactor = 0.1;
    immutable double minLeafLineLength;
    
    this(Cluster* clusters) {       
        this.clusters = clusters; 
        this.leftLim = -0.5;
        int nLeaves = 0;
        
        foreach(leaf; *clusters) {
            _names ~= leaf.name;
            nLeaves++;
        }
        
        this.nLeaves = nLeaves;
        
        double minDist = double.infinity, maxDist = -double.infinity;
        void doDist(Cluster* node) {
            if(node.isLeaf) return;
            minDist = min(node.distance, minDist);
            maxDist = max(node.distance, maxDist);
            
            doDist(node.left);
            doDist(node.right);
        }
        
        doDist(clusters);
        
        this.rightLim = nLeaves - 0.5;
        this.upperLim = maxDist;
    
        auto dummyForm = new Figure();
    
        // Find lower limit by adding in the space taken by the leaf names.
        this.minLeafLineLength = leafLineFactor * (maxDist - minDist);
        this.lowerLim = minDist - minLeafLineLength;
    }
    
protected:
    override void drawPlot(
        Figure form,
        double leftMargin,
        double topMargin,
        double plotWidth,
        double plotHeight
    ) {
        mixin(toPixels);
        int curLeafIndex = 0;
        
        // TODO:  Make these settings customizable.
        auto pen = form.getPen(getColor(0, 0, 0), 1);
        
        double drawLeaf(double parentDistance) {
            immutable parentDistPixels = toPixelsY(parentDistance);
            immutable lineBottomPixels = toPixelsY(lowerLim);
            immutable xPixels = toPixelsX(curLeafIndex);

            form.drawClippedLine(pen,
                PlotPoint(xPixels, lineBottomPixels),
                PlotPoint(xPixels, parentDistPixels)
            );
            
            immutable double ret = curLeafIndex;
            curLeafIndex++;
            return ret;
        }          
        
        // Returns:  The center of the just-drawn sub-dendrogram.
        double drawImpl(Cluster* node) {
            double leftCenter, rightCenter;
            if(node.left.isLeaf) {
                leftCenter = drawLeaf(node.distance);
            } else {
                leftCenter = drawImpl(node.left);
            }
            
            if(node.right.isLeaf) {
                rightCenter = drawLeaf(node.distance);
            } else {
                rightCenter = drawImpl(node.right);
            }
            
            immutable leftTop = (node.left.isLeaf) ? (node.distance) :
                node.left.distance;
            immutable rightTop = (node.right.isLeaf) ? (node.distance) :
                node.right.distance;                
            
            immutable leftCenterPixels = toPixelsX(leftCenter);
            immutable rightCenterPixels = toPixelsX(rightCenter);
            immutable leftTopPixels = toPixelsY(leftTop);
            immutable rightTopPixels = toPixelsY(rightTop);
            immutable distPixels = toPixelsY(node.distance);
            
            form.drawClippedLine(pen, 
                PlotPoint(leftCenterPixels, leftTopPixels),
                PlotPoint(leftCenterPixels, distPixels), 
            );
            
            form.drawClippedLine(pen,
                PlotPoint(leftCenterPixels, distPixels),
                PlotPoint(rightCenterPixels, distPixels)
            );
            
            form.drawClippedLine(pen,
                PlotPoint(rightCenterPixels, distPixels),
                PlotPoint(rightCenterPixels, rightTopPixels)
            );
            
            return 0.5 * leftCenter + 0.5 * rightCenter;
        }
        
        drawImpl(clusters);
    }
    
    protected override void drawLegendSymbol(FigureBase fig, PlotRect where) {
        auto pen = fig.getPen(getColor(0, 0, 0), 1);
        scope(exit) doneWith(pen);
        drawLineLegend(pen, fig, where);
    }
            
public:
    /** 
    Create a Dendrogram from hierarchically clustered data.
    */
    static Dendrogram opCall(Cluster* clusters) {
        return new Dendrogram(clusters);
    }
    
    /**
    Create a Dendrogram figure with rotated text leaf labels on the X-axis
    tick labels and "Distance" as the Y-axis label.
    */
    override Figure toLabeledFigure() @property {
        auto ret = this.toFigure;

        ret.xTickLabels(iota(nLeaves), _names);    
        ret.yLabel("Distance");
        ret.rotatedXTick(true);
        
        return ret;
    }
    
    /**
    The names of the leaf nodes, in order from left to right on the plot.
    */
    final const(string)[] leafNames() const pure nothrow {
        return _names;
    }        
}     

private:

void addFudgeFactors(P)(P plot) {
    with(plot) {
        // Add fudge factors to bounds to make points not appear off the chart.
        immutable horizFudge = (upperLim - lowerLim) * 0.03;
        immutable verticalFudge = (rightLim - leftLim) * 0.03;
        leftLim -= verticalFudge;
        rightLim += verticalFudge;
        upperLim += horizFudge;
        lowerLim -= horizFudge;
    }
}

void drawLineLegend(Pen pen, FigureBase fig, PlotRect where) {
    auto mid = where.y + where.height / 2;
    fig.drawLine(pen,
        PlotPoint(where.x, mid),
        PlotPoint(where.x + where.width, mid));
}

void drawFillLegend(Color color, FigureBase fig, PlotRect where) {
    auto brush = fig.getBrush(color);
    scope(exit) doneWith(brush);
    fig.fillRectangle(brush, where.x, where.y, where.width, where.height);
}

void drawTextLegend(char symbol, Color color, FigureBase fig, PlotRect where) {
    // Draw this small even if the actual symbol is huge so it fits in the
    // allotted space.
    auto font = getFont(plot2kill.util.defaultFont,
        10 + Figure.fontSizeAdjust);

    // Center location.
    string writeThis = [cast(immutable) symbol];
    immutable meas = fig.measureText(writeThis, font);
    immutable stdMeas = fig.measureText("A", font);

    where.y += (where.height - meas.height) / 2 - (stdMeas.height - meas.height);
    where.height = meas.height;

    scope(exit) doneWith(font);
    fig.drawText(writeThis, font, color, where, TextAlignment.Center);
}