analysis/app: add trend graphing
This is still not feature complete, so the page is not linked from
anywhere. It supports graphing benchmark results that match any
perf query. One line is shown per matching benchmark name. Demo
available at:
https://quentin.perf.golang.org/trend?q=upload:20170127.25
Big TODO tasks:
- improve display of multiple benchmarks (multiple graphs?
click-to-select? additional user controls?)
- correlate benchmarks across uploads (currently all uploads are
assumed to be from identical machines)
- optimize/cache queries (6 months of a single benchmark loads in 12s;
6 months of all benchmarks times out at >60s)
Change-Id: I7d1f6073c4837ef63205a10dd4e18085fe7d00ac
Reviewed-on: https://go-review.googlesource.com/36612
Run-TryBot: Quentin Smith <quentin@golang.org>
Reviewed-by: Russ Cox <rsc@golang.org>
TryBot-Result: Gobot Gobot <gobot@golang.org>
diff --git a/analysis/app/app.go b/analysis/app/app.go
index 5ab9459..c4d3029 100644
--- a/analysis/app/app.go
+++ b/analysis/app/app.go
@@ -23,6 +23,7 @@
mux.HandleFunc("/", a.index)
mux.HandleFunc("/search", a.search)
mux.HandleFunc("/compare", a.compare)
+ mux.HandleFunc("/trend", a.trend)
}
// search handles /search.
diff --git a/analysis/app/kza.go b/analysis/app/kza.go
new file mode 100644
index 0000000..9ad6fa6
--- /dev/null
+++ b/analysis/app/kza.go
@@ -0,0 +1,125 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package app
+
+import "math"
+
+// TODO: This all assumes that data is sampled at a regular interval
+// and there are no missing values. It could be generalized to accept
+// missing values (perhaps represented by NaN), or generalized much
+// further by accepting (t, x) pairs and a vector of times at which to
+// evaluate the filter (and an arbitrary window size). I would have to
+// figure out how that affects the difference array in KZA.
+
+// TODO: These can generate a lot of garbage. Perhaps the caller
+// should pass in the target slice? Or these should just overwrite the
+// input array and leave it to the caller to copy if necessary?
+
+// MovingAverage performs a moving average (MA) filter of xs with
+// window size m. m must be a positive odd integer.
+//
+// Note that this is filter is often described in terms of the half
+// length of the window (m-1)/2.
+func MovingAverage(xs []float64, m int) []float64 {
+ if m <= 0 || m%2 != 1 {
+ panic("m must be a positive, odd integer")
+ }
+ ys := make([]float64, len(xs))
+ sum, n := 0.0, 0
+ for l, i, r := -m, -(m-1)/2, 0; i < len(ys); l, i, r = l+1, i+1, r+1 {
+ if l >= 0 {
+ sum -= xs[l]
+ n--
+ }
+ if r < len(xs) {
+ sum += xs[r]
+ n++
+ }
+ if i >= 0 {
+ ys[i] = sum / float64(n)
+ }
+ }
+ return ys
+}
+
+// KolmogorovZurbenko performs a Kolmogorov-Zurbenko (KZ) filter of xs
+// with window size m and k iterations. m must be a positive odd
+// integer. k must be positive.
+func KolmogorovZurbenko(xs []float64, m, k int) []float64 {
+ // k is typically small, and MA is quite efficient, so just do
+ // the iterated moving average rather than bothering to
+ // compute the binomial coefficient kernel.
+ for i := 0; i < k; i++ {
+ // TODO: Generate less garbage.
+ xs = MovingAverage(xs, m)
+ }
+ return xs
+}
+
+// AdaptiveKolmogorovZurbenko performs an adaptive Kolmogorov-Zurbenko
+// (KZA) filter of xs using an initial window size m and k iterations.
+// m must be a positive odd integer. k must be positive.
+//
+// See Zurbenko, et al. 1996: Detecting discontinuities in time series
+// of upper air data: Demonstration of an adaptive filter technique.
+// Journal of Climate, 9, 3548–3560.
+func AdaptiveKolmogorovZurbenko(xs []float64, m, k int) []float64 {
+ // Perform initial KZ filter.
+ z := KolmogorovZurbenko(xs, m, k)
+
+ // Compute differenced values.
+ q := (m - 1) / 2
+ d := make([]float64, len(z)+1)
+ maxD := 0.0
+ for i := q; i < len(z)-q; i++ {
+ d[i] = math.Abs(z[i+q] - z[i-q])
+ if d[i] > maxD {
+ maxD = d[i]
+ }
+ }
+
+ if maxD == 0 {
+ // xs is constant, so no amount of filtering will do
+ // anything. Avoid dividing 0/0 below.
+ return xs
+ }
+
+ // Compute adaptive filter.
+ ys := make([]float64, len(xs))
+ for t := range ys {
+ dPrime := d[t+1] - d[t]
+ f := 1 - d[t]/maxD
+
+ qt := q
+ if dPrime <= 0 {
+ // Zurbenko doesn't specify what to do with
+ // the fractional part of qt and qh, so we
+ // interpret this as summing all points of xs
+ // between qt and qh.
+ qt = int(math.Ceil(float64(q) * f))
+ }
+ if t-qt < 0 {
+ qt = t
+ }
+
+ qh := q
+ if dPrime >= 0 {
+ qh = int(math.Floor(float64(q) * f))
+ }
+ if t+qh >= len(xs) {
+ qh = len(xs) - t - 1
+ }
+
+ sum := 0.0
+ for i := t - qt; i <= t+qh; i++ {
+ sum += xs[i]
+ }
+ // Zurbenko divides by qh+qt, but this undercounts the
+ // number of terms in the sum by 1.
+ ys[t] = sum / float64(qh+qt+1)
+ }
+
+ return ys
+}
diff --git a/analysis/app/kza_test.go b/analysis/app/kza_test.go
new file mode 100644
index 0000000..df65223
--- /dev/null
+++ b/analysis/app/kza_test.go
@@ -0,0 +1,54 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package app
+
+import (
+ "math/rand"
+ "testing"
+)
+
+// Aeq returns true if expect and got are equal to 8 significant
+// figures (1 part in 100 million).
+func Aeq(expect, got float64) bool {
+ if expect < 0 && got < 0 {
+ expect, got = -expect, -got
+ }
+ return expect*0.99999999 <= got && got*0.99999999 <= expect
+}
+
+func TestMovingAverage(t *testing.T) {
+ // Test MovingAverage against the obvious (but slow)
+ // implementation.
+ xs := make([]float64, 100)
+ for iter := 0; iter < 10; iter++ {
+ for i := range xs {
+ xs[i] = rand.Float64()
+ }
+ m := 1 + 2*rand.Intn(100)
+ ys1, ys2 := MovingAverage(xs, m), slowMovingAverage(xs, m)
+
+ // TODO: Use stuff from mathtest.
+ for i, y1 := range ys1 {
+ if !Aeq(y1, ys2[i]) {
+ t.Fatalf("want %v, got %v", ys2, ys1)
+ }
+ }
+ }
+}
+
+func slowMovingAverage(xs []float64, m int) []float64 {
+ ys := make([]float64, len(xs))
+ for i := range ys {
+ psum, n := 0.0, 0
+ for j := i - (m-1)/2; j <= i+(m-1)/2; j++ {
+ if 0 <= j && j < len(xs) {
+ psum += xs[j]
+ n++
+ }
+ }
+ ys[i] = psum / float64(n)
+ }
+ return ys
+}
diff --git a/analysis/app/trend.go b/analysis/app/trend.go
new file mode 100644
index 0000000..cf4837f
--- /dev/null
+++ b/analysis/app/trend.go
@@ -0,0 +1,410 @@
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Loosely based on github.com/aclements/go-misc/benchplot
+
+package app
+
+import (
+ "bytes"
+ "encoding/json"
+ "fmt"
+ "html/template"
+ "io/ioutil"
+ "math"
+ "net/http"
+ "sort"
+ "strconv"
+ "strings"
+
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/ggstat"
+ "github.com/aclements/go-gg/table"
+ "golang.org/x/net/context"
+ "golang.org/x/perf/storage"
+)
+
+// trend handles /trend.
+// With no query, it prints the list of recent uploads containg a "trend" key.
+// With a query, it shows a graph of the matching benchmark results.
+func (a *App) trend(w http.ResponseWriter, r *http.Request) {
+ ctx := requestContext(r)
+
+ if err := r.ParseForm(); err != nil {
+ http.Error(w, err.Error(), 500)
+ return
+ }
+
+ q := r.Form.Get("q")
+
+ tmpl, err := ioutil.ReadFile("template/trend.html")
+ if err != nil {
+ http.Error(w, err.Error(), 500)
+ return
+ }
+
+ t, err := template.New("main").Parse(string(tmpl))
+ if err != nil {
+ http.Error(w, err.Error(), 500)
+ return
+ }
+
+ data := a.trendQuery(ctx, q)
+
+ w.Header().Set("Content-Type", "text/html; charset=utf-8")
+ if err := t.Execute(w, data); err != nil {
+ http.Error(w, err.Error(), 500)
+ return
+ }
+}
+
+// trendData is the struct passed to the trend.html template.
+type trendData struct {
+ Q string
+ Error string
+ TrendUploads []storage.UploadInfo
+ PlotData template.JS
+}
+
+// trendData computes the values for the template and returns a trendData for display.
+func (a *App) trendQuery(ctx context.Context, q string) *trendData {
+ d := &trendData{Q: q}
+ if q == "" {
+ ul := a.StorageClient.ListUploads(`trend>`, []string{"by", "upload-time", "trend"}, 16)
+ defer ul.Close()
+ for ul.Next() {
+ d.TrendUploads = append(d.TrendUploads, ul.Info())
+ }
+ if err := ul.Err(); err != nil {
+ errorf(ctx, "failed to fetch recent trend uploads: %v", err)
+ }
+ return d
+ }
+
+ // TODO(quentin): Chunk query based on matching upload IDs.
+ res := a.StorageClient.Query(q)
+ defer res.Close()
+ t, resultCols := queryToTable(res)
+ for _, col := range []string{"commit", "commit-time", "branch"} {
+ if t.Column(col) == nil {
+ d.Error = fmt.Sprintf("results missing %q label", col)
+ return d
+ }
+ }
+ data := plot(t, resultCols)
+
+ // TODO(quentin): Give the user control over across vs. plotting in separate graphs, instead of only showing one graph with ns/op for each benchmark.
+
+ // Pivot all of the benchmarks into columns of a single table.
+ ar := &aggResults{
+ Across: "name",
+ Values: []string{"filtered normalized mean result", "normalized mean result", "normalized min result", "normalized max result"},
+ }
+ data = ggstat.Agg("commit", "branch", "commit-index")(ar.agg).F(data)
+
+ tables := data.Tables()
+ infof(ctx, "tables: %v", tables)
+ columns := []column{{Name: "commit-index"}}
+ for _, prefix := range ar.Prefixes {
+ if len(ar.Prefixes) == 1 {
+ columns = append(columns,
+ column{Name: prefix + "/normalized mean result"},
+ column{Name: prefix + "/normalized min result", Role: "interval"},
+ column{Name: prefix + "/normalized max result", Role: "interval"},
+ )
+ }
+ columns = append(columns,
+ column{Name: prefix + "/filtered normalized mean result"},
+ )
+ }
+ d.PlotData = tableToJS(data.Table(tables[0]), columns)
+ return d
+}
+
+// queryToTable converts the result of a Query into a Table for later processing.
+// Each label is placed in a column named after the key.
+// Each metric is placed in a separate result column named after the unit.
+func queryToTable(q *storage.Query) (t *table.Table, resultCols []string) {
+ var names []string
+ labels := make(map[string][]string)
+ results := make(map[string][]float64)
+ i := 0
+ for q.Next() {
+ res := q.Result()
+ // TODO(quentin): Handle multiple results with the same name but different NameLabels.
+ names = append(names, res.NameLabels["name"])
+ for k := range res.Labels {
+ if labels[k] == nil {
+ labels[k] = make([]string, i)
+ }
+ }
+ for k := range labels {
+ labels[k] = append(labels[k], res.Labels[k])
+ }
+ f := strings.Fields(res.Content)
+ metrics := make(map[string]float64)
+ for j := 2; j+2 <= len(f); j += 2 {
+ val, err := strconv.ParseFloat(f[j], 64)
+ if err != nil {
+ continue
+ }
+ unit := f[j+1]
+ if results[unit] == nil {
+ results[unit] = make([]float64, i)
+ }
+ metrics[unit] = val
+ }
+ for k := range results {
+ results[k] = append(results[k], metrics[k])
+ }
+ i++
+ }
+
+ tab := new(table.Builder).Add("name", names)
+
+ for k, v := range labels {
+ tab.Add(k, v)
+ }
+ for k, v := range results {
+ tab.Add(k, v)
+ resultCols = append(resultCols, k)
+ }
+
+ sort.Strings(resultCols)
+
+ return tab.Done(), resultCols
+}
+
+// plot takes raw benchmark data in t and produces a Grouping object containing filtered, normalized metric results for a graph.
+// t must contain the string columns "commit", "commit-time", "branch". resultCols specifies the names of float64 columns containing metric results.
+// The returned grouping has columns "commit", "commit-time", "commit-index", "branch", "metric", "normalized min result", "normalized max result", "normalized mean result", "filtered normalized mean result".
+// This is roughly the algorithm from github.com/aclements/go-misc/benchplot
+func plot(t table.Grouping, resultCols []string) table.Grouping {
+ nrows := len(table.GroupBy(t, "name").Tables())
+
+ // Turn ordered commit-time into a "commit-index" column.
+ t = table.SortBy(t, "commit-time", "commit")
+ t = commitIndex{}.F(t)
+
+ // Unpivot all of the metrics into one column.
+ t = table.Unpivot(t, "metric", "result", resultCols...)
+
+ // TODO(quentin): Let user choose which metric(s) to keep.
+ t = table.FilterEq(t, "metric", "ns/op")
+
+ // Average each result at each commit (but keep columns names
+ // the same to keep things easier to read).
+ t = ggstat.Agg("commit", "name", "metric", "branch", "commit-index")(ggstat.AggMean("result"), ggstat.AggMin("result"), ggstat.AggMax("result")).F(t)
+ y := "mean result"
+
+ // Normalize to earliest commit on master. It's important to
+ // do this before the geomean if there are commits missing.
+ // Unfortunately, that also means we have to *temporarily*
+ // group by name and metric, since the geomean needs to be
+ // done on a different grouping.
+ t = table.GroupBy(t, "name", "metric")
+ t = ggstat.Normalize{X: "branch", By: firstMasterIndex, Cols: []string{"mean result", "max result", "min result"}, DenomCols: []string{"mean result", "mean result", "mean result"}}.F(t)
+ y = "normalized " + y
+ for _, col := range []string{"mean result", "max result", "min result"} {
+ t = table.Remove(t, col)
+ }
+ t = table.Ungroup(table.Ungroup(t))
+
+ // Compute geomean for each metric at each commit if there's
+ // more than one benchmark.
+ if len(table.GroupBy(t, "name").Tables()) > 1 {
+ gt := removeNaNs(t, y)
+ gt = ggstat.Agg("commit", "metric", "branch", "commit-index")(ggstat.AggGeoMean(y), ggstat.AggMin("normalized min result"), ggstat.AggMax("normalized max result")).F(gt)
+ gt = table.MapTables(gt, func(_ table.GroupID, t *table.Table) *table.Table {
+ return table.NewBuilder(t).AddConst("name", " geomean").Done()
+ })
+ gt = table.Rename(gt, "geomean "+y, y)
+ gt = table.Rename(gt, "min normalized min result", "normalized min result")
+ gt = table.Rename(gt, "max normalized max result", "normalized max result")
+ t = table.Concat(t, gt)
+ nrows++
+ }
+
+ // Filter the data to reduce noise.
+ t = table.GroupBy(t, "name", "metric")
+ t = kza{y, 15, 3}.F(t)
+ y = "filtered " + y
+ t = table.Ungroup(table.Ungroup(t))
+
+ return t
+}
+
+// aggResults pivots the table, taking the columns in Values and making a new column for each distinct value in Across.
+// aggResults("in", []string{"value1", "value2"} will reshape a table like
+// in value1 value2
+// one 1 2
+// two 3 4
+// and will turn in into a table like
+// one/value1 one/value2 two/value1 two/value2
+// 1 2 3 4
+// across columns must be []string, and value columns must be []float64.
+type aggResults struct {
+ // Across is the name of the column whose values are the column prefix.
+ Across string
+ // Values is the name of the columns to split.
+ Values []string
+ // Prefixes is filled in after calling agg with the name of each prefix that was found.
+ Prefixes []string
+}
+
+// agg implements ggstat.Aggregator and allows using a with ggstat.Agg.
+func (a *aggResults) agg(input table.Grouping, output *table.Builder) {
+ var prefixes []string
+ rows := len(input.Tables())
+ columns := make(map[string][]float64)
+ for i, gid := range input.Tables() {
+ var vs [][]float64
+ for _, col := range a.Values {
+ vs = append(vs, input.Table(gid).MustColumn(col).([]float64))
+ }
+ as := input.Table(gid).MustColumn(a.Across).([]string)
+ for j, prefix := range as {
+ for k, col := range a.Values {
+ key := prefix + "/" + col
+ if columns[key] == nil {
+ if k == 0 {
+ // First time we've seen this prefix, track it.
+ prefixes = append(prefixes, prefix)
+ }
+ columns[key] = make([]float64, rows)
+ for i := range columns[key] {
+ columns[key][i] = math.NaN()
+ }
+ }
+ columns[key][i] = vs[k][j]
+ }
+ }
+ }
+ sort.Strings(prefixes)
+ a.Prefixes = prefixes
+ for _, prefix := range prefixes {
+ for _, col := range a.Values {
+ key := prefix + "/" + col
+ output.Add(key, columns[key])
+ }
+ }
+}
+
+// firstMasterIndex returns the index of the first commit on master.
+// This is used to find the value to normalize against.
+func firstMasterIndex(bs []string) int {
+ return slice.Index(bs, "master")
+}
+
+// commitIndex is a gg.Stat that adds a column called "commit-index" sequentially counting unique values of the column "commit".
+type commitIndex struct{}
+
+func (commitIndex) F(g table.Grouping) table.Grouping {
+ return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+ idxs := make([]int, t.Len())
+ last, idx := "", -1
+ for i, hash := range t.MustColumn("commit").([]string) {
+ if hash != last {
+ idx++
+ last = hash
+ }
+ idxs[i] = idx
+ }
+ t = table.NewBuilder(t).Add("commit-index", idxs).Done()
+
+ return t
+ })
+}
+
+// removeNaNs returns a new Grouping with rows containg NaN in col removed.
+func removeNaNs(g table.Grouping, col string) table.Grouping {
+ return table.Filter(g, func(result float64) bool {
+ return !math.IsNaN(result)
+ }, col)
+}
+
+// kza implements adaptive Kolmogorov-Zurbenko filtering on the data in X.
+type kza struct {
+ X string
+ M, K int
+}
+
+func (k kza) F(g table.Grouping) table.Grouping {
+ return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+ var xs []float64
+ slice.Convert(&xs, t.MustColumn(k.X))
+ nxs := AdaptiveKolmogorovZurbenko(xs, k.M, k.K)
+ return table.NewBuilder(t).Add("filtered "+k.X, nxs).Done()
+ })
+}
+
+// column represents a column in a google.visualization.DataTable
+type column struct {
+ Name string `json:"id"`
+ Role string `json:"role,omitempty"`
+ // These fields are filled in by tableToJS if unspecified.
+ Type string `json:"type"`
+ Label string `json:"label"`
+}
+
+// tableToJS converts a Table to a javascript literal which can be passed to "new google.visualization.DataTable".
+func tableToJS(t *table.Table, columns []column) template.JS {
+ var out bytes.Buffer
+ fmt.Fprint(&out, "{cols: [")
+ var slices []table.Slice
+ for i, c := range columns {
+ if i > 0 {
+ fmt.Fprint(&out, ",\n")
+ }
+ col := t.Column(c.Name)
+ slices = append(slices, col)
+ if c.Type == "" {
+ switch col.(type) {
+ case []string:
+ c.Type = "string"
+ case []int, []float64:
+ c.Type = "number"
+ }
+ }
+ if c.Label == "" {
+ c.Label = c.Name
+ }
+ data, err := json.Marshal(c)
+ if err != nil {
+ panic(err)
+ }
+ out.Write(data)
+ }
+ fmt.Fprint(&out, "],\nrows: [")
+ for i := 0; i < t.Len(); i++ {
+ if i > 0 {
+ fmt.Fprint(&out, ",\n")
+ }
+ fmt.Fprint(&out, "{c:[")
+ for j := range columns {
+ if j > 0 {
+ fmt.Fprint(&out, ", ")
+ }
+ fmt.Fprint(&out, "{v: ")
+ var value []byte
+ var err error
+ switch column := slices[j].(type) {
+ case []string:
+ value, err = json.Marshal(column[i])
+ case []int:
+ value, err = json.Marshal(column[i])
+ case []float64:
+ value, err = json.Marshal(column[i])
+ }
+ if err != nil {
+ panic(err)
+ }
+ out.Write(value)
+ fmt.Fprint(&out, "}")
+ }
+ fmt.Fprint(&out, "]}")
+ }
+ fmt.Fprint(&out, "]}")
+ return template.JS(out.String())
+}
diff --git a/analysis/app/trend_test.go b/analysis/app/trend_test.go
new file mode 100644
index 0000000..8bf2157
--- /dev/null
+++ b/analysis/app/trend_test.go
@@ -0,0 +1,29 @@
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package app
+
+import (
+ "testing"
+
+ "github.com/aclements/go-gg/table"
+ "golang.org/x/perf/internal/diff"
+)
+
+func TestTableToJS(t *testing.T) {
+ in := table.TableFromStrings(
+ []string{"text", "num"},
+ [][]string{
+ {"hello", "15.1"},
+ {"world", "20"},
+ }, true)
+ have := tableToJS(in, []column{{Name: "text"}, {Name: "num"}})
+ want := `{cols: [{"id":"text","type":"string","label":"text"},
+{"id":"num","type":"number","label":"num"}],
+rows: [{c:[{v: "hello"}, {v: 15.1}]},
+{c:[{v: "world"}, {v: 20}]}]}`
+ if d := diff.Diff(string(have), want); d != "" {
+ t.Errorf("tableToJS returned wrong JS (- have, + want):\n%s", d)
+ }
+}
diff --git a/analysis/appengine/template/trend.html b/analysis/appengine/template/trend.html
new file mode 100644
index 0000000..cd1d8eb
--- /dev/null
+++ b/analysis/appengine/template/trend.html
@@ -0,0 +1,73 @@
+<!DOCTYPE html>
+<html>
+ <head>
+ <title>Performance Result Comparison</title>
+ <script type="text/javascript" src="https://www.gstatic.com/charts/loader.js"></script>
+ <style type="text/css">
+#header h1 {
+ display: inline;
+}
+#search {
+ padding: 1em .5em;
+ width: 100%;
+}
+input[type="text"] {
+ font-size: 100%;
+}
+#results {
+ border-top: 1px solid black;
+}
+ </style>
+ </head>
+ <body>
+ <div id="header">
+ <h1>Go Performance Dashboard</h1>
+ <a href="/">about</a>
+ </div>
+ <div id="search">
+ <form action="/trend">
+ <input type="text" name="q" value="{{.Q}}" size="120">
+ <input type="submit" value="Search">
+ </form>
+ </div>
+ <div id="results">
+ {{if not .Q}}
+ <h2>Recent Uploads</h2>
+ <table>
+ <tr><th>Upload ID</th><th>trend</th></tr>
+ {{range .TrendUploads}}
+ <tr><td><a href="/trend?q=upload:{{.UploadID}}">{{.UploadID}}</a></td><td>{{.LabelValues.trend}}</td></tr>
+ {{end}}
+ </table>
+ {{else}}
+ {{with .Error}}
+ <p>{{.}}</p>
+ {{else}}
+ <div id="chart" style="height: 600px"></div>
+ <script type="text/javascript">
+ google.charts.load('current', {'packages':['corechart']});
+ google.charts.setOnLoadCallback(draw);
+ function draw() {
+ var dt = new google.visualization.DataTable({{.PlotData}});
+ var options = {
+ title: 'Benchmark Trend',
+ hAxis: {
+ title: 'commit index',
+ },
+ vAxis: {
+ title: 'normalized ns/op',
+ },
+ explorer: {
+ actions: ['dragToZoom', 'rightClickToReset'],
+ maxZoomIn: 0.05,
+ },
+ };
+ var chart = new google.visualization.LineChart(document.getElementById('chart'));
+ chart.draw(dt, options);
+ }
+ </script>
+ {{end}}
+ {{end}}
+ </div>
+ </body>
+</html>
diff --git a/vendor/github.com/aclements/go-gg/LICENSE b/vendor/github.com/aclements/go-gg/LICENSE
new file mode 100644
index 0000000..a5389da
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/LICENSE
@@ -0,0 +1,28 @@
+Copyright (c) 2016 The Go Authors. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+ * Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+ * Neither the name of Google Inc. nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
diff --git a/vendor/github.com/aclements/go-gg/generic/doc.go b/vendor/github.com/aclements/go-gg/generic/doc.go
new file mode 100644
index 0000000..0df3e8a
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/doc.go
@@ -0,0 +1,6 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package generic provides type-generic functions.
+package generic
diff --git a/vendor/github.com/aclements/go-gg/generic/error.go b/vendor/github.com/aclements/go-gg/generic/error.go
new file mode 100644
index 0000000..a04e5e5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/error.go
@@ -0,0 +1,21 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package generic
+
+import "reflect"
+
+type TypeError struct {
+ Type1, Type2 reflect.Type
+ Extra string
+}
+
+func (e TypeError) Error() string {
+ msg := e.Type1.String()
+ if e.Type2 != nil {
+ msg += " and " + e.Type2.String()
+ }
+ msg += " " + e.Extra
+ return msg
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/order.go b/vendor/github.com/aclements/go-gg/generic/order.go
new file mode 100644
index 0000000..551373d
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/order.go
@@ -0,0 +1,96 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package generic
+
+import "reflect"
+
+// CanOrder returns whether the values a and b are orderable according
+// to the Go language specification.
+func CanOrder(a, b interface{}) bool {
+ ak, bk := reflect.ValueOf(a).Kind(), reflect.ValueOf(b).Kind()
+ if ak != bk {
+ return false
+ }
+ return CanOrderR(ak)
+}
+
+var orderable = map[reflect.Kind]bool{
+ reflect.Int: true,
+ reflect.Int8: true,
+ reflect.Int16: true,
+ reflect.Int32: true,
+ reflect.Int64: true,
+ reflect.Uint: true,
+ reflect.Uintptr: true,
+ reflect.Uint8: true,
+ reflect.Uint16: true,
+ reflect.Uint32: true,
+ reflect.Uint64: true,
+ reflect.Float32: true,
+ reflect.Float64: true,
+ reflect.String: true,
+}
+
+// CanOrderR returns whether two values of kind k are orderable
+// according to the Go language specification.
+func CanOrderR(k reflect.Kind) bool {
+ return orderable[k]
+}
+
+// Order returns the order of values a and b: -1 if a < b, 0 if a ==
+// b, 1 if a > b. The results are undefined if either a or b is NaN.
+//
+// Order panics if a and b are not orderable according to the Go
+// language specification.
+func Order(a, b interface{}) int {
+ return OrderR(reflect.ValueOf(a), reflect.ValueOf(b))
+}
+
+// OrderR is equivalent to Order, but operates on reflect.Values.
+func OrderR(a, b reflect.Value) int {
+ if a.Kind() != b.Kind() {
+ panic(&TypeError{a.Type(), b.Type(), "are not orderable because they are different kinds"})
+ }
+
+ switch a.Kind() {
+ case reflect.Float32, reflect.Float64:
+ a, b := a.Float(), b.Float()
+ if a < b {
+ return -1
+ } else if a > b {
+ return 1
+ }
+ return 0
+
+ case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+ a, b := a.Int(), b.Int()
+ if a < b {
+ return -1
+ } else if a > b {
+ return 1
+ }
+ return 0
+
+ case reflect.Uint, reflect.Uintptr, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64:
+ a, b := a.Uint(), b.Uint()
+ if a < b {
+ return -1
+ } else if a > b {
+ return 1
+ }
+ return 0
+
+ case reflect.String:
+ a, b := a.String(), b.String()
+ if a < b {
+ return -1
+ } else if a > b {
+ return 1
+ }
+ return 0
+ }
+
+ panic(&TypeError{a.Type(), nil, "is not orderable"})
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/concat.go b/vendor/github.com/aclements/go-gg/generic/slice/concat.go
new file mode 100644
index 0000000..e538790
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/concat.go
@@ -0,0 +1,41 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Concat returns the concatenation of all of ss. The types of all of
+// the arguments must be identical or Concat will panic with a
+// *generic.TypeError. The returned slice will have the same type as the
+// inputs. If there are 0 arguments, Concat returns nil. Concat does
+// not modify any of the input slices.
+func Concat(ss ...T) T {
+ if len(ss) == 0 {
+ return nil
+ }
+
+ rvs := make([]reflect.Value, len(ss))
+ total := 0
+ var typ reflect.Type
+ for i, s := range ss {
+ rvs[i] = reflectSlice(s)
+ total += rvs[i].Len()
+ if i == 0 {
+ typ = rvs[i].Type()
+ } else if rvs[i].Type() != typ {
+ panic(&generic.TypeError{typ, rvs[i].Type(), "have different types"})
+ }
+ }
+
+ out := reflect.MakeSlice(typ, 0, total)
+ for _, rv := range rvs {
+ out = reflect.AppendSlice(out, rv)
+ }
+ return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/convert.go b/vendor/github.com/aclements/go-gg/generic/slice/convert.go
new file mode 100644
index 0000000..c3d46a3
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/convert.go
@@ -0,0 +1,55 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Convert converts each element in from and assigns it to *to. to
+// must be a pointer to a slice. Convert slices or extends *to to
+// len(from) and then assigns to[i] = T(from[i]) where T is the type
+// of *to's elements. If from and *to have the same element type, it
+// simply assigns *to = from.
+func Convert(to interface{}, from T) {
+ fv := reflectSlice(from)
+ tv := reflect.ValueOf(to)
+ if tv.Kind() != reflect.Ptr {
+ panic(&generic.TypeError{tv.Type(), nil, "is not a *[]T"})
+ }
+ tst := tv.Type().Elem()
+ if tst.Kind() != reflect.Slice {
+ panic(&generic.TypeError{tv.Type(), nil, "is not a *[]T"})
+ }
+
+ if fv.Type().AssignableTo(tst) {
+ tv.Elem().Set(fv)
+ return
+ }
+
+ eltt := tst.Elem()
+ if !fv.Type().Elem().ConvertibleTo(eltt) {
+ panic(&generic.TypeError{fv.Type(), tst, "cannot be converted"})
+ }
+
+ switch to := to.(type) {
+ case *[]float64:
+ // This is extremely common.
+ *to = (*to)[:0]
+ for i, len := 0, fv.Len(); i < len; i++ {
+ *to = append(*to, fv.Index(i).Convert(eltt).Float())
+ }
+
+ default:
+ tsv := tv.Elem()
+ tsv.SetLen(0)
+ for i, len := 0, fv.Len(); i < len; i++ {
+ tsv = reflect.Append(tsv, fv.Index(i).Convert(eltt))
+ }
+ tv.Elem().Set(tsv)
+ }
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/cycle.go b/vendor/github.com/aclements/go-gg/generic/slice/cycle.go
new file mode 100644
index 0000000..a94e5e2
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/cycle.go
@@ -0,0 +1,45 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import "reflect"
+
+// Cycle constructs a slice of length length by repeatedly
+// concatenating s to itself. If len(s) >= length, it returns
+// s[:length]. Otherwise, it allocates a new slice. If len(s) == 0 and
+// length != 0, Cycle panics.
+func Cycle(s T, length int) T {
+ rv := reflectSlice(s)
+ if rv.Len() >= length {
+ return rv.Slice(0, length).Interface()
+ }
+
+ if rv.Len() == 0 {
+ panic("empty slice")
+ }
+
+ // Allocate a new slice of the appropriate length.
+ out := reflect.MakeSlice(rv.Type(), length, length)
+
+ // Copy elements to out.
+ for pos := 0; pos < length; {
+ pos += reflect.Copy(out.Slice(pos, length), rv)
+ }
+
+ return out.Interface()
+}
+
+// Repeat returns a slice consisting of length copies of v.
+func Repeat(v interface{}, length int) T {
+ if length < 0 {
+ length = 0
+ }
+ rv := reflect.ValueOf(v)
+ out := reflect.MakeSlice(reflect.SliceOf(rv.Type()), length, length)
+ for i := 0; i < length; i++ {
+ out.Index(i).Set(rv)
+ }
+ return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/doc.go b/vendor/github.com/aclements/go-gg/generic/slice/doc.go
new file mode 100644
index 0000000..2ebb862
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/doc.go
@@ -0,0 +1,6 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package slice provides generic slice functions.
+package slice
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/find.go b/vendor/github.com/aclements/go-gg/generic/slice/find.go
new file mode 100644
index 0000000..09c7d53
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/find.go
@@ -0,0 +1,51 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Index returns the index of the first instance of val in s, or -1 if
+// val is not present in s. val's type must be s's element type.
+func Index(s T, val interface{}) int {
+ rs := reflectSlice(s)
+ if vt := reflect.TypeOf(val); rs.Type().Elem() != vt {
+ // TODO: Better "<seq> is not a sequence of <val>".
+ panic(&generic.TypeError{rs.Type(), vt, "cannot find"})
+ }
+
+ for i, l := 0, rs.Len(); i < l; i++ {
+ if rs.Index(i).Interface() == val {
+ return i
+ }
+ }
+ return -1
+}
+
+// LastIndex returns the index of the last instance of val in s, or -1
+// if val is not present in s. val's type must be s's element type.
+func LastIndex(s T, val interface{}) int {
+ rs := reflectSlice(s)
+ if vt := reflect.TypeOf(val); rs.Type().Elem() != vt {
+ // TODO: Better "<seq> is not a sequence of <val>".
+ panic(&generic.TypeError{rs.Type(), vt, "cannot find"})
+ }
+
+ for i := rs.Len() - 1; i >= 0; i-- {
+ if rs.Index(i).Interface() == val {
+ return i
+ }
+ }
+ return -1
+}
+
+// Contains reports whether val is within s. val's type must be s's
+// element type.
+func Contains(s T, val interface{}) bool {
+ return Index(s, val) >= 0
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/index.go b/vendor/github.com/aclements/go-gg/generic/slice/index.go
new file mode 100644
index 0000000..780b4ac
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/index.go
@@ -0,0 +1,82 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Select returns a slice w such that w[i] = v[indexes[i]].
+func Select(v T, indexes []int) T {
+ switch v := v.(type) {
+ case []int:
+ res := make([]int, len(indexes))
+ for i, x := range indexes {
+ res[i] = v[x]
+ }
+ return res
+
+ case []float64:
+ res := make([]float64, len(indexes))
+ for i, x := range indexes {
+ res[i] = v[x]
+ }
+ return res
+
+ case []string:
+ res := make([]string, len(indexes))
+ for i, x := range indexes {
+ res[i] = v[x]
+ }
+ return res
+ }
+
+ rv := reflectSlice(v)
+ res := reflect.MakeSlice(rv.Type(), len(indexes), len(indexes))
+ for i, x := range indexes {
+ res.Index(i).Set(rv.Index(x))
+ }
+ return res.Interface()
+}
+
+// SelectInto assigns out[i] = in[indexes[i]]. in and out must have
+// the same types and len(out) must be >= len(indexes). If in and out
+// overlap, the results are undefined.
+func SelectInto(out, in T, indexes []int) {
+ // TODO: Maybe they should only have to be assignable?
+ if it, ot := reflect.TypeOf(in), reflect.TypeOf(out); it != ot {
+ panic(&generic.TypeError{it, ot, "must be the same type"})
+ }
+
+ switch in := in.(type) {
+ case []int:
+ out := out.([]int)
+ for i, x := range indexes {
+ out[i] = in[x]
+ }
+ return
+
+ case []float64:
+ out := out.([]float64)
+ for i, x := range indexes {
+ out[i] = in[x]
+ }
+ return
+
+ case []string:
+ out := out.([]string)
+ for i, x := range indexes {
+ out[i] = in[x]
+ }
+ return
+ }
+
+ inv, outv := reflectSlice(in), reflectSlice(out)
+ for i, x := range indexes {
+ outv.Index(i).Set(inv.Index(x))
+ }
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/min.go b/vendor/github.com/aclements/go-gg/generic/slice/min.go
new file mode 100644
index 0000000..7cda3fc
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/min.go
@@ -0,0 +1,99 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+ "sort"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Min returns the minimum value in v. v must either implement
+// sort.Interface or its elements must be orderable. Min panics if v
+// is empty.
+func Min(v T) interface{} {
+ x, _ := minmax(v, -1, true)
+ return x.Interface()
+}
+
+// ArgMin returns the index of the minimum value in v. If there are
+// multiple indexes equal to the minimum value, ArgMin returns the
+// lowest of them. v must be a slice whose elements are orderable, or
+// must implement sort.Interface. ArgMin panics if v is empty.
+func ArgMin(v interface{}) int {
+ _, i := minmax(v, -1, false)
+ return i
+}
+
+// Max returns the maximum value in v. v must either implement
+// sort.Interface or its elements must be orderable. Max panics if v
+// is empty.
+func Max(v T) interface{} {
+ x, _ := minmax(v, 1, true)
+ return x.Interface()
+}
+
+// ArgMax returns the index of the maximum value in v. If there are
+// multiple indexes equal to the maximum value, ArgMax returns the
+// lowest of them. v must be a slice whose elements are orderable, or
+// must implement sort.Interface. ArgMax panics if v is empty.
+func ArgMax(v interface{}) int {
+ _, i := minmax(v, 1, false)
+ return i
+}
+
+func minmax(v interface{}, keep int, val bool) (reflect.Value, int) {
+ switch v := v.(type) {
+ case sort.Interface:
+ if v.Len() == 0 {
+ if keep < 0 {
+ panic("zero-length sequence has no minimum")
+ } else {
+ panic("zero-length sequence has no maximum")
+ }
+ }
+ maxi := 0
+ if keep < 0 {
+ for i, len := 0, v.Len(); i < len; i++ {
+ if v.Less(i, maxi) {
+ maxi = i
+ }
+ }
+ } else {
+ for i, len := 0, v.Len(); i < len; i++ {
+ if v.Less(maxi, i) {
+ maxi = i
+ }
+ }
+ }
+
+ if !val {
+ return reflect.Value{}, maxi
+ }
+
+ rv := reflectSlice(v)
+ return rv.Index(maxi), maxi
+ }
+
+ rv := reflectSlice(v)
+ if !generic.CanOrderR(rv.Type().Elem().Kind()) {
+ panic(&generic.TypeError{rv.Type().Elem(), nil, "is not orderable"})
+ }
+ if rv.Len() == 0 {
+ if keep < 0 {
+ panic("zero-length slice has no minimum")
+ } else {
+ panic("zero-length slice has no maximum")
+ }
+ }
+ max, maxi := rv.Index(0), 0
+ for i, len := 1, rv.Len(); i < len; i++ {
+ if elt := rv.Index(i); generic.OrderR(elt, max) == keep {
+ max, maxi = elt, i
+ }
+ }
+ return max, maxi
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/nub.go b/vendor/github.com/aclements/go-gg/generic/slice/nub.go
new file mode 100644
index 0000000..0c5f69c
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/nub.go
@@ -0,0 +1,52 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import "reflect"
+
+var trueVal = reflect.ValueOf(true)
+
+// Nub returns v with duplicates removed. It keeps the first instance
+// of each distinct value and preserves their order.
+func Nub(v T) T {
+ rv := reflectSlice(v)
+ indexes := make([]int, 0)
+ set := reflect.MakeMap(reflect.MapOf(rv.Type().Elem(), trueVal.Type()))
+ for i, l := 0, rv.Len(); i < l; i++ {
+ x := rv.Index(i)
+ if set.MapIndex(x).IsValid() {
+ continue
+ }
+ set.SetMapIndex(x, trueVal)
+ indexes = append(indexes, i)
+ }
+ return Select(v, indexes)
+}
+
+// NubAppend is equivalent to appending all of the slices in vs and
+// then calling Nub on the result, but more efficient.
+func NubAppend(vs ...T) T {
+ if len(vs) == 0 {
+ return nil
+ }
+
+ rv := reflectSlice(vs[0])
+ set := reflect.MakeMap(reflect.MapOf(rv.Type().Elem(), trueVal.Type()))
+ out := reflect.MakeSlice(rv.Type(), 0, 0)
+
+ for _, v := range vs {
+ rv := reflectSlice(v)
+ for i, l := 0, rv.Len(); i < l; i++ {
+ x := rv.Index(i)
+ if set.MapIndex(x).IsValid() {
+ continue
+ }
+ set.SetMapIndex(x, trueVal)
+ out = reflect.Append(out, x)
+ }
+ }
+
+ return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/seq.go b/vendor/github.com/aclements/go-gg/generic/slice/seq.go
new file mode 100644
index 0000000..47812e4
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/seq.go
@@ -0,0 +1,28 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// T is a Go slice value of type []U.
+//
+// This is primarily for documentation. There is no way to statically
+// enforce this in Go; however, functions that expect a slice will
+// panic with a *generic.TypeError if passed a non-slice value.
+type T interface{}
+
+// reflectSlice checks that s is a slice and returns its
+// reflect.Value. It panics with a *generic.TypeError if s is not a slice.
+func reflectSlice(s T) reflect.Value {
+ rv := reflect.ValueOf(s)
+ if rv.Kind() != reflect.Slice {
+ panic(&generic.TypeError{rv.Type(), nil, "is not a slice"})
+ }
+ return rv
+}
diff --git a/vendor/github.com/aclements/go-gg/generic/slice/sort.go b/vendor/github.com/aclements/go-gg/generic/slice/sort.go
new file mode 100644
index 0000000..e5ef8b6
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/generic/slice/sort.go
@@ -0,0 +1,137 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package slice
+
+import (
+ "reflect"
+ "sort"
+ "time"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// CanSort returns whether the value v can be sorted.
+func CanSort(v interface{}) bool {
+ switch v.(type) {
+ case sort.Interface, []time.Time:
+ return true
+ }
+ return generic.CanOrderR(reflect.TypeOf(v).Elem().Kind())
+}
+
+// Sort sorts v in increasing order. v must implement sort.Interface,
+// be a slice whose elements are orderable, or be a []time.Time.
+func Sort(v interface{}) {
+ sort.Sort(Sorter(v))
+}
+
+// Sorter returns a sort.Interface for sorting v. v must implement
+// sort.Interface, be a slice whose elements are orderable, or be a
+// []time.Time.
+func Sorter(v interface{}) sort.Interface {
+ switch v := v.(type) {
+ case []int:
+ return sort.IntSlice(v)
+ case []float64:
+ return sort.Float64Slice(v)
+ case []string:
+ return sort.StringSlice(v)
+ case []time.Time:
+ return sortTimeSlice(v)
+ case sort.Interface:
+ return v
+ }
+
+ rv := reflectSlice(v)
+ switch rv.Type().Elem().Kind() {
+ case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+ return sortIntSlice{rv}
+ case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr:
+ return sortUintSlice{rv}
+ case reflect.Float32, reflect.Float64:
+ return sortFloatSlice{rv}
+ case reflect.String:
+ return sortStringSlice{rv}
+ }
+ panic(&generic.TypeError{rv.Type().Elem(), nil, "is not orderable"})
+}
+
+type sortIntSlice struct {
+ reflect.Value
+}
+
+func (s sortIntSlice) Len() int {
+ return s.Value.Len()
+}
+
+func (s sortIntSlice) Less(i, j int) bool {
+ return s.Index(i).Int() < s.Index(j).Int()
+}
+
+func (s sortIntSlice) Swap(i, j int) {
+ a, b := s.Index(i).Int(), s.Index(j).Int()
+ s.Index(i).SetInt(b)
+ s.Index(j).SetInt(a)
+}
+
+type sortUintSlice struct {
+ reflect.Value
+}
+
+func (s sortUintSlice) Len() int {
+ return s.Value.Len()
+}
+
+func (s sortUintSlice) Less(i, j int) bool {
+ return s.Index(i).Uint() < s.Index(j).Uint()
+}
+
+func (s sortUintSlice) Swap(i, j int) {
+ a, b := s.Index(i).Uint(), s.Index(j).Uint()
+ s.Index(i).SetUint(b)
+ s.Index(j).SetUint(a)
+}
+
+type sortFloatSlice struct {
+ reflect.Value
+}
+
+func (s sortFloatSlice) Len() int {
+ return s.Value.Len()
+}
+
+func (s sortFloatSlice) Less(i, j int) bool {
+ return s.Index(i).Float() < s.Index(j).Float()
+}
+
+func (s sortFloatSlice) Swap(i, j int) {
+ a, b := s.Index(i).Float(), s.Index(j).Float()
+ s.Index(i).SetFloat(b)
+ s.Index(j).SetFloat(a)
+}
+
+type sortStringSlice struct {
+ reflect.Value
+}
+
+func (s sortStringSlice) Len() int {
+ return s.Value.Len()
+}
+
+func (s sortStringSlice) Less(i, j int) bool {
+ return s.Index(i).String() < s.Index(j).String()
+}
+
+func (s sortStringSlice) Swap(i, j int) {
+ a, b := s.Index(i).String(), s.Index(j).String()
+ s.Index(i).SetString(b)
+ s.Index(j).SetString(a)
+}
+
+type sortTimeSlice []time.Time
+
+func (s sortTimeSlice) Len() int { return len(s) }
+func (s sortTimeSlice) Less(i, j int) bool { return s[i].Before(s[j]) }
+func (s sortTimeSlice) Swap(i, j int) { s[i], s[j] = s[j], s[i] }
diff --git a/vendor/github.com/aclements/go-gg/ggstat/agg.go b/vendor/github.com/aclements/go-gg/ggstat/agg.go
new file mode 100644
index 0000000..072d621
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/agg.go
@@ -0,0 +1,297 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "fmt"
+ "reflect"
+
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/stats"
+ "github.com/aclements/go-moremath/vec"
+)
+
+// TODO: AggFirst, AggTukey. StdDev?
+
+// Agg constructs an Aggregate transform from a grouping column and a
+// set of Aggregators.
+//
+// TODO: Does this belong in ggstat? The specific aggregator functions
+// probably do, but the concept could go in package table.
+func Agg(xs ...string) func(aggs ...Aggregator) Aggregate {
+ return func(aggs ...Aggregator) Aggregate {
+ return Aggregate{xs, aggs}
+ }
+}
+
+// Aggregate computes aggregate functions of a table grouped by
+// distinct values of a column or set of columns.
+//
+// Aggregate first groups the table by the Xs columns. Each of these
+// groups produces a single row in the output table, where the unique
+// value of each of the Xs columns appears in the output row, along
+// with constant columns from the input, as well as any columns that
+// have a unique value within every group (they're "effectively"
+// constant). Additional columns in the output row are produced by
+// applying the Aggregator functions to the group.
+type Aggregate struct {
+ // Xs is the list column names to group values by before
+ // computing aggregate functions.
+ Xs []string
+
+ // Aggregators is the set of Aggregator functions to apply to
+ // each group of values.
+ Aggregators []Aggregator
+}
+
+// An Aggregator is a function that aggregates each group of input
+// into one row and adds it to output. It may be based on multiple
+// columns from input and may add multiple columns to output.
+type Aggregator func(input table.Grouping, output *table.Builder)
+
+func (s Aggregate) F(g table.Grouping) table.Grouping {
+ isConst := make([]bool, len(g.Columns()))
+ for i := range isConst {
+ isConst[i] = true
+ }
+
+ subgroups := map[table.GroupID]table.Grouping{}
+ for _, gid := range g.Tables() {
+ g := table.GroupBy(g.Table(gid), s.Xs...)
+ subgroups[gid] = g
+
+ for i, col := range g.Columns() {
+ if !isConst[i] {
+ continue
+ }
+ // Can this column be promoted to constant?
+ for _, gid2 := range g.Tables() {
+ t := g.Table(gid2)
+ isConst[i] = isConst[i] && checkConst(t, col)
+ }
+ }
+ }
+
+ return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+ g := table.GroupBy(t, s.Xs...)
+ var nt table.Builder
+
+ // Construct X columns.
+ rows := len(g.Tables())
+ for colidx, xcol := range s.Xs {
+ xs := reflect.MakeSlice(table.ColType(t, xcol), rows, rows)
+ for i, gid := range g.Tables() {
+ for j := 0; j < len(s.Xs)-colidx-1; j++ {
+ gid = gid.Parent()
+ }
+ xs.Index(i).Set(reflect.ValueOf(gid.Label()))
+ }
+
+ nt.Add(xcol, xs.Interface())
+ }
+
+ // Apply Aggregators.
+ for _, agg := range s.Aggregators {
+ agg(g, &nt)
+ }
+
+ // Keep constant and effectively constant columns.
+ for i := range isConst {
+ col := t.Columns()[i]
+ if !isConst[i] || nt.Has(col) {
+ continue
+ }
+ if cv, ok := t.Const(col); ok {
+ nt.AddConst(col, cv)
+ continue
+ }
+
+ ncol := reflect.MakeSlice(table.ColType(t, col), len(g.Tables()), len(g.Tables()))
+ for i, gid := range g.Tables() {
+ v := reflect.ValueOf(g.Table(gid).Column(col))
+ ncol.Index(i).Set(v.Index(0))
+ }
+ nt.Add(col, ncol.Interface())
+ }
+ return nt.Done()
+ })
+}
+
+func checkConst(t *table.Table, col string) bool {
+ if _, ok := t.Const(col); ok {
+ return true
+ }
+ v := reflect.ValueOf(t.Column(col))
+ if v.Len() <= 1 {
+ return true
+ }
+ if !v.Type().Elem().Comparable() {
+ return false
+ }
+ elem := v.Index(0).Interface()
+ for i, l := 1, v.Len(); i < l; i++ {
+ if elem != v.Index(i).Interface() {
+ return false
+ }
+ }
+ return true
+}
+
+// AggCount returns an aggregate function that computes the number of
+// rows in each group. The resulting column will be named label, or
+// "count" if label is "".
+func AggCount(label string) Aggregator {
+ if label == "" {
+ label = "count"
+ }
+
+ return func(input table.Grouping, b *table.Builder) {
+ counts := make([]int, 0, len(input.Tables()))
+ for _, gid := range input.Tables() {
+ counts = append(counts, input.Table(gid).Len())
+ }
+ b.Add(label, counts)
+ }
+}
+
+// AggMean returns an aggregate function that computes the mean of
+// each of cols. The resulting columns will be named "mean <col>" and
+// will have the same type as <col>.
+func AggMean(cols ...string) Aggregator {
+ return aggFn(stats.Mean, "mean ", cols...)
+}
+
+// AggGeoMean returns an aggregate function that computes the
+// geometric mean of each of cols. The resulting columns will be named
+// "geomean <col>" and will have the same type as <col>.
+func AggGeoMean(cols ...string) Aggregator {
+ return aggFn(stats.GeoMean, "geomean ", cols...)
+}
+
+// AggMin returns an aggregate function that computes the minimum of
+// each of cols. The resulting columns will be named "min <col>" and
+// will have the same type as <col>.
+func AggMin(cols ...string) Aggregator {
+ min := func(xs []float64) float64 {
+ x, _ := stats.Bounds(xs)
+ return x
+ }
+ return aggFn(min, "min ", cols...)
+}
+
+// AggMax returns an aggregate function that computes the maximum of
+// each of cols. The resulting columns will be named "max <col>" and
+// will have the same type as <col>.
+func AggMax(cols ...string) Aggregator {
+ max := func(xs []float64) float64 {
+ _, x := stats.Bounds(xs)
+ return x
+ }
+ return aggFn(max, "max ", cols...)
+}
+
+// AggSum returns an aggregate function that computes the sum of each
+// of cols. The resulting columns will be named "sum <col>" and will
+// have the same type as <col>.
+func AggSum(cols ...string) Aggregator {
+ return aggFn(vec.Sum, "sum ", cols...)
+}
+
+// AggQuantile returns an aggregate function that computes a quantile
+// of each of cols. quantile has a range of [0,1]. The resulting
+// columns will be named "<prefix> <col>" and will have the same type
+// as <col>.
+func AggQuantile(prefix string, quantile float64, cols ...string) Aggregator {
+ // "prefix" could be autogenerated (e.g. fmt.Sprintf("p%g ",
+ // quantile * 100)), but then the caller would need to do the
+ // same fmt.Sprintf to compute the column name they had just
+ // created. Perhaps Aggregator should provide a way to find
+ // the generated column names.
+ return aggFn(func(data []float64) float64 {
+ return stats.Sample{Xs: data}.Quantile(quantile)
+ }, prefix+" ", cols...)
+}
+
+func aggFn(f func([]float64) float64, prefix string, cols ...string) Aggregator {
+ ocols := make([]string, len(cols))
+ for i, col := range cols {
+ ocols[i] = prefix + col
+ }
+
+ return func(input table.Grouping, b *table.Builder) {
+ for coli, col := range cols {
+ means := make([]float64, 0, len(input.Tables()))
+
+ var xs []float64
+ var ct reflect.Type
+ for i, gid := range input.Tables() {
+ v := input.Table(gid).MustColumn(col)
+ if i == 0 {
+ ct = reflect.TypeOf(v)
+ }
+ slice.Convert(&xs, v)
+ means = append(means, f(xs))
+ }
+
+ if ct == float64SliceType {
+ b.Add(ocols[coli], means)
+ } else {
+ // Convert means back to the type of col.
+ outptr := reflect.New(ct)
+ slice.Convert(outptr.Interface(), means)
+ b.Add(ocols[coli], outptr.Elem().Interface())
+ }
+ }
+ }
+}
+
+// AggUnique returns an aggregate function retains the unique value of
+// each of cols within each aggregate group, or panics if some group
+// contains more than one value for one of these columns.
+//
+// Note that Aggregate will automatically retain columns that happen
+// to be unique. AggUnique can be used to enforce at aggregation time
+// that certain columns *must* be unique (and get a nice error if they
+// are not).
+func AggUnique(cols ...string) Aggregator {
+ return func(input table.Grouping, b *table.Builder) {
+ if len(cols) == 0 {
+ return
+ }
+ if len(input.Tables()) == 0 {
+ panic(fmt.Sprintf("unknown column: %q", cols[0]))
+ }
+
+ for _, col := range cols {
+ ctype := table.ColType(input, col)
+ rows := len(input.Tables())
+ vs := reflect.MakeSlice(ctype, rows, rows)
+ for i, gid := range input.Tables() {
+ // Get values in this column.
+ xs := reflect.ValueOf(input.Table(gid).MustColumn(col))
+
+ // Check for uniqueness.
+ if xs.Len() == 0 {
+ panic(fmt.Sprintf("cannot AggUnique empty column %q", col))
+ }
+ uniquev := xs.Index(0)
+ unique := uniquev.Interface()
+ for i, len := 1, xs.Len(); i < len; i++ {
+ other := xs.Index(i).Interface()
+ if unique != other {
+ panic(fmt.Sprintf("column %q is not unique; contains at least %v and %v", col, unique, other))
+ }
+ }
+
+ // Store unique value.
+ vs.Index(i).Set(uniquev)
+ }
+
+ // Add unique values slice to output table.
+ b.Add(col, vs.Interface())
+ }
+ }
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/bin.go b/vendor/github.com/aclements/go-gg/ggstat/bin.go
new file mode 100644
index 0000000..2602004
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/bin.go
@@ -0,0 +1,196 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "math"
+ "reflect"
+ "sort"
+
+ "github.com/aclements/go-gg/generic"
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/vec"
+)
+
+// XXX If this is just based on the number of bins, it can come up
+// with really ugly boundary numbers. If the bin width is specified,
+// then you could also specify the left edge and bins will be placed
+// at [align+width*N, align+width*(N+1)]. ggplot2 also lets you
+// specify the center alignment.
+//
+// XXX In Matlab and NumPy, bins are open on the right *except* for
+// the last bin, which is closed on both.
+//
+// XXX Number of bins/bin width/specify boundaries, same bins across
+// all groups/separate for each group/based on shared scales (don't
+// have that information here), relative or absolute histogram (Matlab
+// has lots more).
+//
+// XXX Scale transform.
+//
+// The result of Bin has two columns in addition to constant columns from the input:
+//
+// - Column X is the left edge of the bin.
+//
+// - Column W is the sum of the rows' weights, or column "count" is
+// the number of rows in the bin.
+type Bin struct {
+ // X is the name of the column to use for samples.
+ X string
+
+ // W is the optional name of the column to use for sample
+ // weights. It may be "" to weight each sample as 1.
+ W string
+
+ // Width controls how wide each bin should be. If not provided
+ // or 0, a width will be chosen to produce 30 bins. If X is an
+ // integer column, this width will be treated as an integer as
+ // well.
+ Width float64
+
+ // Center controls the center point of each bin. To center on
+ // integers, for example, you could use {Width: 1, Center:
+ // 0}.
+ // XXX What does center mean for integers? Should an unspecified center yield an autochosen one, or 0?
+ //Center float64
+
+ // Breaks is the set of break points to use as boundaries
+ // between bins. The interval of each bin is [Breaks[i],
+ // Breaks[i+1]). Data points before the first break are
+ // dropped. If provided, Width and Center are ignored.
+ Breaks table.Slice
+
+ // SplitGroups indicates that each group in the table should
+ // have separate bounds based on the data in that group alone.
+ // The default, false, indicates that the binning function
+ // should use the bounds of all of the data combined. This
+ // makes it easier to compare bins across groups.
+ SplitGroups bool
+}
+
+func (b Bin) F(g table.Grouping) table.Grouping {
+ breaks := reflect.ValueOf(b.Breaks)
+ agg := AggCount("count")
+ if b.W != "" {
+ agg = aggFn(vec.Sum, "", b.W)
+ }
+ if !breaks.IsValid() && !b.SplitGroups {
+ breaks = b.computeBreaks(g)
+ }
+ // Change b.X to the start of the bin.
+ g = table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+ breaks := breaks
+ if !breaks.IsValid() {
+ breaks = b.computeBreaks(t)
+ }
+ nbreaks := breaks.Len()
+
+ in := reflect.ValueOf(t.MustColumn(b.X))
+ nin := in.Len()
+
+ out := reflect.MakeSlice(breaks.Type(), nin, nin)
+ var found []int
+ for i := 0; i < nin; i++ {
+ elt := in.Index(i)
+ bin := sort.Search(nbreaks, func(j int) bool {
+ return generic.OrderR(elt, breaks.Index(j)) < 0
+ })
+ // 0 means the row doesn't fit on the front
+ // XXX Allow configuring the first and last bin as infinite or not.
+ bin = bin - 1
+ if bin >= 0 {
+ found = append(found, i)
+ out.Index(i).Set(breaks.Index(bin))
+ }
+ }
+ var nt table.Builder
+ for _, col := range t.Columns() {
+ if col == b.X {
+ nt.Add(col, slice.Select(out.Interface(), found))
+ } else if c, ok := t.Const(col); ok {
+ nt.AddConst(col, c)
+ } else {
+ nt.Add(col, slice.Select(t.Column(col), found))
+ }
+ }
+ return nt.Done()
+ })
+ // Group by the found bin
+ return Agg(b.X)(agg).F(g)
+}
+
+func (b Bin) computeBreaks(g table.Grouping) reflect.Value {
+ var cols []slice.T
+ for _, gid := range g.Tables() {
+ cols = append(cols, g.Table(gid).MustColumn(b.X))
+ }
+ data := slice.Concat(cols...)
+
+ min := slice.Min(data)
+ max := slice.Max(data)
+
+ rv := reflect.ValueOf(min)
+ switch rv.Type().Kind() {
+ case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64:
+ min, max := rv.Int(), reflect.ValueOf(max).Int()
+ width := int64(b.Width)
+ if width == 0 {
+ width = (max - min) / 30
+ if width < 1 {
+ width = 1
+ }
+ }
+ // XXX: This assumes boundaries should be aligned with
+ // 0. We should support explicit Center or Boundary
+ // requests.
+ min -= (min % width)
+ var breaks []int64
+ for i := min; i < max; i += width {
+ breaks = append(breaks, i)
+ }
+ outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+ slice.Convert(outs.Interface(), breaks)
+ return outs.Elem()
+ case reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr:
+ min, max := rv.Uint(), reflect.ValueOf(max).Uint()
+ width := uint64(b.Width)
+ if width == 0 {
+ width = (max - min) / 30
+ if width < 1 {
+ width = 1
+ }
+ }
+ min -= (min % width)
+ var breaks []uint64
+ for i := min; i < max; i += width {
+ breaks = append(breaks, i)
+ }
+ outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+ slice.Convert(outs.Interface(), breaks)
+ return outs.Elem()
+ case reflect.Float32, reflect.Float64:
+ min, max := rv.Float(), reflect.ValueOf(max).Float()
+ width := b.Width
+ if width == 0 {
+ width = (max - min) / 30
+ if width == 0 {
+ width = 1
+ }
+ }
+ min -= math.Mod(min, width)
+ var breaks []float64
+ for i := min; i < max; i += width {
+ breaks = append(breaks, i)
+ }
+ outs := reflect.New(reflect.ValueOf(cols[0]).Type())
+ slice.Convert(outs.Interface(), breaks)
+ return outs.Elem()
+ default:
+ panic("can't compute breaks for unknown type")
+ }
+}
+
+// TODO: Count for categorical data.
diff --git a/vendor/github.com/aclements/go-gg/ggstat/common.go b/vendor/github.com/aclements/go-gg/ggstat/common.go
new file mode 100644
index 0000000..db0d7a5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/common.go
@@ -0,0 +1,10 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import "reflect"
+
+var float64Type = reflect.TypeOf(float64(0))
+var float64SliceType = reflect.TypeOf([]float64(nil))
diff --git a/vendor/github.com/aclements/go-gg/ggstat/density.go b/vendor/github.com/aclements/go-gg/ggstat/density.go
new file mode 100644
index 0000000..73b8a60
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/density.go
@@ -0,0 +1,123 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/stats"
+ "github.com/aclements/go-moremath/vec"
+)
+
+// TODO: Default to first (and second) column for X (and Y)?
+
+// Density constructs a probability density estimate from a set of
+// samples using kernel density estimation.
+//
+// X is the only required field. All other fields have reasonable
+// default zero values.
+//
+// The result of Density has three columns in addition to constant
+// columns from the input:
+//
+// - Column X is the points at which the density estimate is sampled.
+//
+// - Column "probability density" is the density estimate.
+//
+// - Column "cumulative density" is the cumulative density estimate.
+type Density struct {
+ // X is the name of the column to use for samples.
+ X string
+
+ // W is the optional name of the column to use for sample
+ // weights. It may be "" to uniformly weight samples.
+ W string
+
+ // N is the number of points to sample the KDE at. If N is 0,
+ // a reasonable default is used.
+ //
+ // TODO: This is particularly sensitive to the scale
+ // transform.
+ //
+ // TODO: Base the default on the bandwidth. If the bandwidth
+ // is really narrow, we may need a lot of samples to exceed
+ // the Nyquist rate.
+ N int
+
+ // Domain specifies the domain at which to sample this function.
+ // If Domain is nil, it defaults to DomainData{}.
+ Domain FunctionDomainer
+
+ // Kernel is the kernel to use for the KDE.
+ Kernel stats.KDEKernel
+
+ // Bandwidth is the bandwidth to use for the KDE.
+ //
+ // If this is zero, the bandwidth is computed from the data
+ // using a default bandwidth estimator (currently
+ // stats.BandwidthScott).
+ Bandwidth float64
+
+ // BoundaryMethod is the boundary correction method to use for
+ // the KDE. The default value is BoundaryReflect; however, the
+ // default bounds are effectively +/-inf, which is equivalent
+ // to performing no boundary correction.
+ BoundaryMethod stats.KDEBoundaryMethod
+
+ // [BoundaryMin, BoundaryMax) specify a bounded support for
+ // the KDE. If both are 0 (their default values), they are
+ // treated as +/-inf.
+ //
+ // To specify a half-bounded support, set Min to math.Inf(-1)
+ // or Max to math.Inf(1).
+ BoundaryMin float64
+ BoundaryMax float64
+}
+
+func (d Density) F(g table.Grouping) table.Grouping {
+ kde := stats.KDE{
+ Kernel: d.Kernel,
+ Bandwidth: d.Bandwidth,
+ BoundaryMethod: d.BoundaryMethod,
+ BoundaryMin: d.BoundaryMin,
+ BoundaryMax: d.BoundaryMax,
+ }
+ dname, cname := "probability density", "cumulative density"
+
+ addEmpty := func(out *table.Builder) {
+ out.Add(dname, []float64{})
+ out.Add(cname, []float64{})
+ }
+
+ return Function{
+ X: d.X, N: d.N, Domain: d.Domain,
+ Fn: func(gid table.GroupID, in *table.Table, sampleAt []float64, out *table.Builder) {
+ if len(sampleAt) == 0 {
+ addEmpty(out)
+ return
+ }
+
+ // Get input sample.
+ var sample stats.Sample
+ slice.Convert(&sample.Xs, in.MustColumn(d.X))
+ if d.W != "" {
+ slice.Convert(&sample.Weights, in.MustColumn(d.W))
+ if sample.Weight() == 0 {
+ addEmpty(out)
+ return
+ }
+ }
+
+ // Compute KDE.
+ kde.Sample = sample
+ if d.Bandwidth == 0 {
+ kde.Bandwidth = stats.BandwidthScott(sample)
+ }
+
+ out.Add(dname, vec.Map(kde.PDF, sampleAt))
+ out.Add(cname, vec.Map(kde.CDF, sampleAt))
+ },
+ }.F(g)
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/domain.go b/vendor/github.com/aclements/go-gg/ggstat/domain.go
new file mode 100644
index 0000000..77bb78e
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/domain.go
@@ -0,0 +1,109 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "math"
+
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/stats"
+)
+
+// A FunctionDomainer computes the domain over which to evaluate a
+// statistical function.
+type FunctionDomainer interface {
+ // FunctionDomain computes the domain of a particular column
+ // within a table. It takes a Grouping and a column in that
+ // Grouping to compute the domain of and returns a function
+ // that computes the domain for a specific group in the
+ // Grouping. This makes it possible for FunctionDomain to
+ // easily compute either Grouping-wide domains, or per-Table
+ // domains.
+ //
+ // The returned domain may be (NaN, NaN) to indicate that
+ // there is no data and the domain is vacuous.
+ FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64)
+}
+
+// DomainFixed is a FunctionDomainer that returns a fixed domain.
+type DomainFixed struct {
+ Min, Max float64
+}
+
+var _ FunctionDomainer = DomainFixed{}
+
+func (r DomainFixed) FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64) {
+ return func(table.GroupID) (min, max float64) {
+ return r.Min, r.Max
+ }
+}
+
+// DomainData is a FunctionDomainer that computes domains based on the
+// bounds of the data.
+type DomainData struct {
+ // Widen expands the domain by Widen times the span of the
+ // data.
+ //
+ // A value of 1.0 means to use exactly the bounds of the data.
+ // If Widen is 0, it is treated as 1.1 (that is, widen the
+ // domain by 10%, or 5% on the left and 5% on the right).
+ Widen float64
+
+ // SplitGroups indicates that each group in the table should
+ // have a separate domain based on the data in that group
+ // alone. The default, false, indicates that the domain should
+ // be based on all of the data in the table combined. This
+ // makes it possible to stack functions and easier to compare
+ // them across groups.
+ SplitGroups bool
+}
+
+var _ FunctionDomainer = DomainData{}
+
+const defaultWiden = 1.1
+
+func (r DomainData) FunctionDomain(g table.Grouping, col string) func(gid table.GroupID) (min, max float64) {
+ widen := r.Widen
+ if widen <= 0 {
+ widen = defaultWiden
+ }
+
+ var xs []float64
+ if !r.SplitGroups {
+ // Compute combined bounds.
+ gmin, gmax := math.NaN(), math.NaN()
+ for _, gid := range g.Tables() {
+ t := g.Table(gid)
+ slice.Convert(&xs, t.MustColumn(col))
+ xmin, xmax := stats.Bounds(xs)
+ if xmin < gmin || math.IsNaN(gmin) {
+ gmin = xmin
+ }
+ if xmax > gmax || math.IsNaN(gmax) {
+ gmax = xmax
+ }
+ }
+
+ // Widen bounds.
+ span := gmax - gmin
+ gmin, gmax = gmin-span*(widen-1)/2, gmax+span*(widen-1)/2
+
+ return func(table.GroupID) (min, max float64) {
+ return gmin, gmax
+ }
+ }
+
+ return func(gid table.GroupID) (min, max float64) {
+ // Compute bounds.
+ slice.Convert(&xs, g.Table(gid).MustColumn(col))
+ min, max = stats.Bounds(xs)
+
+ // Widen bounds.
+ span := max - min
+ min, max = min-span*(widen-1)/2, max+span*(widen-1)/2
+ return
+ }
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/ecdf.go b/vendor/github.com/aclements/go-gg/ggstat/ecdf.go
new file mode 100644
index 0000000..46364a8
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/ecdf.go
@@ -0,0 +1,138 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/vec"
+)
+
+// ECDF constructs an empirical CDF from a set of samples.
+//
+// X is the only required field. All other fields have reasonable
+// default zero values.
+//
+// The result of ECDF has three columns in addition to constant
+// columns from the input. The names of the columns depend on whether
+// Label is "".
+//
+// - Column X is the points at which the CDF changes (a subset of the
+// samples).
+//
+// - Column "cumulative density" or "cumulative density of <label>" is
+// the cumulative density estimate.
+//
+// - Column "cumulative count" (if W and Label are ""), "cumulative
+// weight" (if W is not "", but Label is "") or "cumulative <label>"
+// (if Label is not "") is the cumulative count or weight of samples.
+// That is, cumulative density times the total weight of the samples.
+type ECDF struct {
+ // X is the name of the column to use for samples.
+ X string
+
+ // W is the optional name of the column to use for sample
+ // weights. It may be "" to uniformly weight samples.
+ W string
+
+ // Label, if not "", gives a label for the samples. It is used
+ // to construct more specific names for the output columns. It
+ // should be a plural noun.
+ Label string
+
+ // Domain specifies the domain of the returned ECDF. If the
+ // domain is wider than the bounds of the data in a group,
+ // ECDF will add a point below the smallest sample and above
+ // the largest sample to make the 0 and 1 levels clear. If
+ // Domain is nil, it defaults to DomainData{}.
+ Domain FunctionDomainer
+}
+
+func (s ECDF) F(g table.Grouping) table.Grouping {
+ // Set defaults.
+ if s.Domain == nil {
+ s.Domain = DomainData{}
+ }
+
+ // Construct output column names.
+ dname, cname := "cumulative density", "cumulative count"
+ if s.Label != "" {
+ dname += " of " + s.Label
+ cname = "cumulative " + s.Label
+ } else if s.W != "" {
+ cname = "cumulative weight"
+ }
+
+ g = table.SortBy(g, s.X)
+ domain := s.Domain.FunctionDomain(g, s.X)
+
+ return table.MapTables(g, func(gid table.GroupID, t *table.Table) *table.Table {
+ // Get input columns.
+ var xs, ws []float64
+ slice.Convert(&xs, t.MustColumn(s.X))
+ if s.W != "" {
+ slice.Convert(&ws, t.MustColumn(s.W))
+ }
+
+ // Ignore empty tables.
+ if len(xs) == 0 {
+ nt := new(table.Builder).Add(s.X, []float64{}).Add(cname, []float64{}).Add(dname, []float64{})
+ preserveConsts(nt, t)
+ return nt.Done()
+ }
+
+ // Get domain.
+ min, max := domain(gid)
+
+ // Create output columns.
+ xo, do, co := make([]float64, 0), make([]float64, 0), make([]float64, 0)
+ if min < xs[0] {
+ // Extend to the left.
+ xo = append(xo, min)
+ do = append(do, 0)
+ co = append(co, 0)
+ }
+
+ // Compute total weight.
+ var total float64
+ if ws == nil {
+ total = float64(t.Len())
+ } else {
+ total = vec.Sum(ws)
+ }
+
+ // Create ECDF.
+ cum := 0.0
+ for i := 0; i < len(xs); {
+ j := i
+ for j < len(xs) && xs[i] == xs[j] {
+ if ws == nil {
+ cum += 1
+ } else {
+ cum += ws[j]
+ }
+ j++
+ }
+
+ xo = append(xo, xs[i])
+ do = append(do, cum/total)
+ co = append(co, cum)
+
+ i = j
+ }
+
+ if xs[len(xs)-1] < max {
+ // Extend to the right.
+ xo = append(xo, max)
+ do = append(do, 1)
+ co = append(co, cum)
+ }
+
+ // Construct results table.
+ nt := new(table.Builder).Add(s.X, xo).Add(dname, do).Add(cname, co)
+ preserveConsts(nt, t)
+ return nt.Done()
+ })
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/fn.go b/vendor/github.com/aclements/go-gg/ggstat/fn.go
new file mode 100644
index 0000000..0a4de33
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/fn.go
@@ -0,0 +1,123 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "math"
+ "reflect"
+
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+ "github.com/aclements/go-moremath/vec"
+)
+
+// Function samples a continuous univariate function at N points in
+// the domain computed by Domain.
+//
+// The result of Function binds column X to the X values at which the
+// function is sampled and retains constant columns from the input.
+// The computed function can add arbitrary columns for its output.
+type Function struct {
+ // X is the name of the column to use for input domain of this
+ // function.
+ X string
+
+ // N is the number of points to sample the function at. If N
+ // is 0, a reasonable default is used.
+ N int
+
+ // Domain specifies the domain of which to sample this function.
+ // If Domain is nil, it defaults to DomainData{}.
+ Domain FunctionDomainer
+
+ // Fn is the continuous univariate function to sample. Fn will
+ // be called with each table in the grouping and the X values
+ // at which it should be sampled. Fn must add its output
+ // columns to out. The output table will already contain the
+ // sample points bound to the X column.
+ Fn func(gid table.GroupID, in *table.Table, sampleAt []float64, out *table.Builder)
+}
+
+const defaultFunctionSamples = 200
+
+func (f Function) F(g table.Grouping) table.Grouping {
+ // Set defaults.
+ if f.N <= 0 {
+ f.N = defaultFunctionSamples
+ }
+ if f.Domain == nil {
+ f.Domain = DomainData{}
+ }
+
+ domain := f.Domain.FunctionDomain(g, f.X)
+ return table.MapTables(g, func(gid table.GroupID, t *table.Table) *table.Table {
+ min, max := domain(gid)
+
+ // Compute sample points. If there's no data, there
+ // are no sample points, but we still have to run the
+ // function to get the right output columns.
+ var ss []float64
+ if math.IsNaN(min) {
+ ss = []float64{}
+ } else {
+ ss = vec.Linspace(min, max, f.N)
+ }
+
+ var nt table.Builder
+ ctype := table.ColType(t, f.X)
+ if ctype == float64Type {
+ // Bind output X column.
+ nt.Add(f.X, ss)
+ } else {
+ // Convert to the column type.
+ vsp := reflect.New(ctype)
+ slice.Convert(vsp.Interface(), ss)
+ vs := vsp.Elem()
+ // This may have produced duplicate values.
+ // Eliminate those.
+ if vs.Len() > 0 {
+ prev, i := vs.Index(0).Interface(), 1
+ for j := 1; j < vs.Len(); j++ {
+ next := vs.Index(j).Interface()
+ if prev == next {
+ // Skip duplicate.
+ continue
+ }
+
+ if i != j {
+ vs.Index(i).Set(vs.Index(j))
+ }
+ i++
+ prev = next
+ }
+ vs.SetLen(i)
+ }
+ // Bind column-typed values to output X.
+ nt.Add(f.X, vs.Interface())
+ // And convert back to []float64 so we can
+ // apply the function.
+ slice.Convert(&ss, vs.Interface())
+ }
+
+ // Apply the function to the sample points.
+ f.Fn(gid, t, ss, &nt)
+
+ preserveConsts(&nt, t)
+ return nt.Done()
+ })
+}
+
+// preserveConsts copies the constant columns from t into nt.
+func preserveConsts(nt *table.Builder, t *table.Table) {
+ for _, col := range t.Columns() {
+ if nt.Has(col) {
+ // Don't overwrite existing columns in nt.
+ continue
+ }
+ if cv, ok := t.Const(col); ok {
+ nt.AddConst(col, cv)
+ }
+ }
+}
diff --git a/vendor/github.com/aclements/go-gg/ggstat/normalize.go b/vendor/github.com/aclements/go-gg/ggstat/normalize.go
new file mode 100644
index 0000000..4644246
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/ggstat/normalize.go
@@ -0,0 +1,175 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package ggstat
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic/slice"
+ "github.com/aclements/go-gg/table"
+)
+
+// Normalize normalizes each group such that some data point is 1.
+//
+// Either X or Index is required (though 0 is a reasonable value of
+// Index).
+//
+// The result of Normalize is the same as the input table, plus
+// additional columns for each normalized column. These columns will
+// be named "normalized <col>" where <col> is the name of the original
+// column and will have type []float64.
+type Normalize struct {
+ // X is the name of the column to use to find the denominator
+ // row. If X is "", Index is used instead.
+ X string
+
+ // Index is the row index of the denominator row if X is ""
+ // (otherwise it is ignored). Index may be negative, in which
+ // case it is added to the number of rows (e.g., -1 is the
+ // last row).
+ Index int
+
+ // By is a function func([]T) int that returns the index of
+ // the denominator row given column X. By may be nil, in which
+ // case it defaults to generic.ArgMin.
+ By interface{}
+
+ // Cols is a slice of the names of columns to normalize
+ // relative to the corresponding DenomCols value in the
+ // denominator row. Cols may be nil, in which case it defaults
+ // to all integral and floating point columns.
+ Cols []string
+
+ // DenomCols is a slice of the names of columns used as the
+ // demoninator. DenomCols may be nil, in which case it
+ // defaults to Cols (i.e. each column will be normalized to
+ // the value from that column in the denominator row.)
+ // Otherwise, DenomCols must be the same length as Cols.
+ DenomCols []string
+}
+
+func (s Normalize) F(g table.Grouping) table.Grouping {
+ // Find the columns to normalize.
+ if s.Cols == nil {
+ cols := []string{}
+ for i, ct := range colTypes(g) {
+ if canNormalize(ct.Elem().Kind()) {
+ cols = append(cols, g.Columns()[i])
+ }
+ }
+ s.Cols = cols
+ }
+ if len(s.Cols) == 0 {
+ return g
+ }
+
+ // Construct new column names.
+ newcols := make([]string, len(s.Cols))
+ for i, col := range s.Cols {
+ newcols[i] = "normalized " + col
+ }
+
+ // Get "by" function.
+ var byv reflect.Value
+ byargs := make([]reflect.Value, 1)
+ if s.By != nil {
+ byv = reflect.ValueOf(s.By)
+ // TODO: Type check byv better.
+ }
+
+ return table.MapTables(g, func(_ table.GroupID, t *table.Table) *table.Table {
+ if t.Len() == 0 {
+ return t
+ }
+
+ // Find the denominator row.
+ var drow int
+ if s.X == "" {
+ drow = s.Index
+ if drow < 0 {
+ drow += t.Len()
+ }
+ } else {
+ xs := t.MustColumn(s.X)
+ if s.By == nil {
+ drow = slice.ArgMin(xs)
+ } else {
+ byargs[0] = reflect.ValueOf(xs)
+ byout := byv.Call(byargs)
+ drow = int(byout[0].Int())
+ }
+ }
+
+ // Normalize columns.
+ newt := table.NewBuilder(t)
+ denomCols := s.DenomCols
+ if denomCols == nil {
+ denomCols = s.Cols
+ }
+ for coli, col := range s.Cols {
+ denom := denomValue(t.MustColumn(denomCols[coli]), drow)
+ out := normalizeTo(t.MustColumn(col), denom)
+ newt.Add(newcols[coli], out)
+ }
+
+ return newt.Done()
+ })
+}
+
+func colTypes(g table.Grouping) []reflect.Type {
+ cts := make([]reflect.Type, len(g.Columns()))
+ for i, col := range g.Columns() {
+ cts[i] = table.ColType(g, col)
+ }
+ return cts
+}
+
+var canNormalizeKinds = map[reflect.Kind]bool{
+ reflect.Float32: true,
+ reflect.Float64: true,
+ reflect.Int: true,
+ reflect.Int8: true,
+ reflect.Int16: true,
+ reflect.Int32: true,
+ reflect.Int64: true,
+ reflect.Uint: true,
+ reflect.Uintptr: true,
+ reflect.Uint8: true,
+ reflect.Uint16: true,
+ reflect.Uint32: true,
+ reflect.Uint64: true,
+}
+
+func canNormalize(k reflect.Kind) bool {
+ return canNormalizeKinds[k]
+}
+
+func denomValue(s interface{}, index int) float64 {
+ switch s := s.(type) {
+ case []float64:
+ return s[index]
+ }
+ return reflect.ValueOf(s).Index(index).Convert(float64Type).Float()
+}
+
+func normalizeTo(s interface{}, denom float64) interface{} {
+ switch s := s.(type) {
+ case []float64:
+ out := make([]float64, len(s))
+ for i, numer := range s {
+ out[i] = numer / denom
+ }
+ return out
+ }
+
+ sv := reflect.ValueOf(s)
+
+ out := reflect.MakeSlice(float64SliceType, sv.Len(), sv.Len())
+ for i, len := 0, sv.Len(); i < len; i++ {
+ numer := sv.Index(i).Convert(float64Type).Float()
+ out.Index(i).SetFloat(numer / denom)
+ }
+ return out.Interface()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/concat.go b/vendor/github.com/aclements/go-gg/table/concat.go
new file mode 100644
index 0000000..14a8f4f
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/concat.go
@@ -0,0 +1,104 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "fmt"
+
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+// Concat returns the concatenation of the rows in each matching group
+// across gs. All Groupings in gs must have the same set of columns
+// (though they need not be in the same order; the column order from
+// gs[0] will be used). The GroupIDs in the returned Grouping will be
+// the union of the GroupIDs in gs.
+func Concat(gs ...Grouping) Grouping {
+ if len(gs) == 0 {
+ return new(Table)
+ }
+
+ // Check that all Groupings have the same set of columns. They
+ // can be in different orders.
+ colSet := map[string]bool{}
+ for _, col := range gs[0].Columns() {
+ colSet[col] = true
+ }
+ for i, g2 := range gs[1:] {
+ diff := len(g2.Columns()) != len(colSet)
+ if !diff {
+ for _, col := range g2.Columns() {
+ if !colSet[col] {
+ diff = true
+ break
+ }
+ }
+ }
+ if diff {
+ panic(fmt.Sprintf("columns in Groupings 0 and %d differ: %q vs %q", i+1, gs[0].Columns(), g2.Columns()))
+ }
+ }
+
+ // Collect group IDs.
+ haveGID := map[GroupID]bool{}
+ gids := []GroupID{}
+ for _, g := range gs {
+ for _, gid := range g.Tables() {
+ if haveGID[gid] {
+ continue
+ }
+ haveGID[gid] = true
+ gids = append(gids, gid)
+ }
+ }
+
+ // Build output groups.
+ var ng GroupingBuilder
+ for _, gid := range gids {
+ // Build output table.
+ var nt Builder
+ var cols []slice.T
+ for _, col := range gs[0].Columns() {
+ // Is it constant?
+ isConst := false
+ var cv interface{}
+ for _, g := range gs {
+ t := g.Table(gid)
+ if t == nil {
+ continue
+ }
+ if cv1, ok := t.Const(col); ok {
+ if !isConst {
+ isConst = true
+ cv = cv1
+ } else if cv != cv1 {
+ isConst = false
+ break
+ }
+ } else {
+ isConst = false
+ break
+ }
+ }
+ if isConst {
+ nt.AddConst(col, cv)
+ continue
+ }
+
+ // Not a constant. Collect slices.
+ for _, g := range gs {
+ t := g.Table(gid)
+ if t == nil {
+ continue
+ }
+ cols = append(cols, t.Column(col))
+ }
+ nt.Add(col, slice.Concat(cols...))
+ cols = cols[:0]
+ }
+ ng.Add(gid, nt.Done())
+ }
+ return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/filter.go b/vendor/github.com/aclements/go-gg/table/filter.go
new file mode 100644
index 0000000..5008138
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/filter.go
@@ -0,0 +1,94 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "fmt"
+ "reflect"
+
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+var boolType = reflect.TypeOf(false)
+
+// Filter filters g to only rows where pred returns true. pred must be
+// a function that returns bool and takes len(cols) arguments where
+// the type of col[i] is assignable to argument i.
+//
+// TODO: Create a faster batch variant where pred takes slices.
+func Filter(g Grouping, pred interface{}, cols ...string) Grouping {
+ // TODO: Use generic.TypeError.
+ predv := reflect.ValueOf(pred)
+ predt := predv.Type()
+ if predt.Kind() != reflect.Func || predt.NumIn() != len(cols) || predt.NumOut() != 1 || predt.Out(0) != boolType {
+ panic("predicate function must be func(col[0], col[1], ...) bool")
+ }
+ if len(cols) == 0 {
+ return g
+ }
+ if len(g.Tables()) == 0 {
+ panic(fmt.Sprintf("unknown column %q", cols[0]))
+ }
+ // Type check arguments.
+ for i, col := range cols {
+ colt := ColType(g, col)
+ if !colt.Elem().AssignableTo(predt.In(i)) {
+ panic(fmt.Sprintf("column %d (type %s) is not assignable to predicate argument %d (type %s)", i, colt.Elem(), i, predt.In(i)))
+ }
+ }
+
+ args := make([]reflect.Value, len(cols))
+ colvs := make([]reflect.Value, len(cols))
+ match := make([]int, 0)
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ // Get columns.
+ for i, col := range cols {
+ colvs[i] = reflect.ValueOf(t.MustColumn(col))
+ }
+
+ // Find the set of row indexes that satisfy pred.
+ match = match[:0]
+ for r, len := 0, t.Len(); r < len; r++ {
+ for c, colv := range colvs {
+ args[c] = colv.Index(r)
+ }
+ if predv.Call(args)[0].Bool() {
+ match = append(match, r)
+ }
+ }
+
+ // Create the new table.
+ if len(match) == t.Len() {
+ return t
+ }
+ var nt Builder
+ for _, col := range t.Columns() {
+ nt.Add(col, slice.Select(t.Column(col), match))
+ }
+ return nt.Done()
+ })
+}
+
+// FilterEq filters g to only rows where the value in col equals val.
+func FilterEq(g Grouping, col string, val interface{}) Grouping {
+ match := make([]int, 0)
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ // Find the set of row indexes that match val.
+ seq := t.MustColumn(col)
+ match = match[:0]
+ rv := reflect.ValueOf(seq)
+ for i, len := 0, rv.Len(); i < len; i++ {
+ if rv.Index(i).Interface() == val {
+ match = append(match, i)
+ }
+ }
+
+ var nt Builder
+ for _, col := range t.Columns() {
+ nt.Add(col, slice.Select(t.Column(col), match))
+ }
+ return nt.Done()
+ })
+}
diff --git a/vendor/github.com/aclements/go-gg/table/group.go b/vendor/github.com/aclements/go-gg/table/group.go
new file mode 100644
index 0000000..afabba2
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/group.go
@@ -0,0 +1,243 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "fmt"
+ "reflect"
+ "strings"
+
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+// GroupID identifies a group. GroupIDs form a tree, rooted at
+// RootGroupID (which is also the zero GroupID).
+type GroupID struct {
+ *groupNode
+}
+
+// RootGroupID is the root of the GroupID tree.
+var RootGroupID = GroupID{}
+
+type groupNode struct {
+ parent GroupID
+ label interface{}
+}
+
+// String returns the path to GroupID g in the form "/l1/l2/l3". If g
+// is RootGroupID, it returns "/". Each level in the group is formed
+// by formatting the label using fmt's "%v" verb. Note that this is
+// purely diagnostic; this string may not uniquely identify g.
+func (g GroupID) String() string {
+ if g == RootGroupID {
+ return "/"
+ }
+ parts := []string{}
+ for p := g; p != RootGroupID; p = p.parent {
+ part := fmt.Sprintf("/%v", p.label)
+ parts = append(parts, part)
+ }
+ for i, j := 0, len(parts)-1; i < j; i, j = i+1, j-1 {
+ parts[i], parts[j] = parts[j], parts[i]
+ }
+ return strings.Join(parts, "")
+}
+
+// Extend returns a new GroupID that is a child of GroupID g. The
+// returned GroupID will not be equal to any existing GroupID (even if
+// label is not unique among g's children). The label is primarily
+// diagnostic; the table package uses it only when printing tables,
+// but callers may store semantic information in group labels.
+func (g GroupID) Extend(label interface{}) GroupID {
+ return GroupID{&groupNode{g, label}}
+}
+
+// Parent returns the parent of g. The parent of RootGroupID is
+// RootGroupID.
+func (g GroupID) Parent() GroupID {
+ if g == RootGroupID {
+ return RootGroupID
+ }
+ return g.parent
+}
+
+// Label returns the label of g.
+func (g GroupID) Label() interface{} {
+ return g.label
+}
+
+// GroupBy sub-divides all groups such that all of the rows in each
+// group have equal values for all of the named columns. The relative
+// order of rows with equal values for the named columns is
+// maintained. Grouped-by columns become constant columns within each
+// group.
+func GroupBy(g Grouping, cols ...string) Grouping {
+ // TODO: This would generate much less garbage if we grouped
+ // all of cols in one pass.
+ //
+ // TODO: This constructs one slice per column per input group,
+ // but it would be even better if it constructed just one
+ // slice per column.
+
+ if len(cols) == 0 {
+ return g
+ }
+
+ var out GroupingBuilder
+ for _, gid := range g.Tables() {
+ t := g.Table(gid)
+
+ if cv, ok := t.Const(cols[0]); ok {
+ // Grouping by a constant is trivial.
+ subgid := gid.Extend(cv)
+ out.Add(subgid, t)
+ continue
+ }
+
+ c := t.MustColumn(cols[0])
+
+ // Create an index on c.
+ type subgroupInfo struct {
+ key interface{}
+ rows []int
+ }
+ subgroups := []subgroupInfo{}
+ keys := make(map[interface{}]int)
+ seq := reflect.ValueOf(c)
+ for i := 0; i < seq.Len(); i++ {
+ x := seq.Index(i).Interface()
+ sg, ok := keys[x]
+ if !ok {
+ sg = len(subgroups)
+ subgroups = append(subgroups, subgroupInfo{x, []int{}})
+ keys[x] = sg
+ }
+ subgroup := &subgroups[sg]
+ subgroup.rows = append(subgroup.rows, i)
+ }
+
+ // Count rows in each subgroup.
+ offsets := make([]int, 1+len(subgroups))
+ for i := range subgroups {
+ offsets[i+1] = offsets[i] + len(subgroups[i].rows)
+ }
+
+ // Split each column.
+ builders := make([]Builder, len(subgroups))
+ for _, name := range t.Columns() {
+ if name == cols[0] {
+ // Promote the group-by column to a
+ // constant.
+ for i := range subgroups {
+ builders[i].AddConst(name, subgroups[i].key)
+ }
+ continue
+ }
+
+ if cv, ok := t.Const(name); ok {
+ // Keep constants constant.
+ for i := range builders {
+ builders[i].AddConst(name, cv)
+ }
+ continue
+ }
+
+ // Create a slice for all of the values.
+ col := t.Column(name)
+ ncol := reflect.MakeSlice(reflect.TypeOf(col), t.Len(), t.Len())
+
+ // Shuffle each subgroup into ncol.
+ for i := range subgroups {
+ subcol := ncol.Slice(offsets[i], offsets[i+1]).Interface()
+ slice.SelectInto(subcol, col, subgroups[i].rows)
+ builders[i].Add(name, subcol)
+ }
+ }
+
+ // Add tables to output Grouping.
+ for i := range builders {
+ subgid := gid.Extend(subgroups[i].key)
+ out.Add(subgid, builders[i].Done())
+ }
+ }
+
+ return GroupBy(out.Done(), cols[1:]...)
+}
+
+// Ungroup concatenates adjacent Tables in g that share a group parent
+// into a Table identified by the parent, undoing the effects of the
+// most recent GroupBy operation.
+func Ungroup(g Grouping) Grouping {
+ groups := g.Tables()
+ if len(groups) == 0 || len(groups) == 1 && groups[0] == RootGroupID {
+ return g
+ }
+
+ var out GroupingBuilder
+ runGid := groups[0].Parent()
+ runTabs := []*Table{}
+ for _, gid := range groups {
+ if gid.Parent() != runGid {
+ // Flush the run.
+ out.Add(runGid, concatRows(runTabs...))
+
+ runGid = gid.Parent()
+ runTabs = runTabs[:0]
+ }
+ runTabs = append(runTabs, g.Table(gid))
+ }
+ // Flush the last run.
+ out.Add(runGid, concatRows(runTabs...))
+
+ return out.Done()
+}
+
+// Flatten concatenates all of the groups in g into a single Table.
+// This is equivalent to repeatedly Ungrouping g.
+func Flatten(g Grouping) *Table {
+ groups := g.Tables()
+ switch len(groups) {
+ case 0:
+ return new(Table)
+
+ case 1:
+ return g.Table(groups[0])
+ }
+
+ tabs := make([]*Table, len(groups))
+ for i, gid := range groups {
+ tabs[i] = g.Table(gid)
+ }
+
+ return concatRows(tabs...)
+}
+
+// concatRows concatenates the rows of tabs into a single Table. All
+// Tables in tabs must all have the same column set.
+func concatRows(tabs ...*Table) *Table {
+ // TODO: Consider making this public. It would have to check
+ // the columns, and we would probably also want a concatCols.
+
+ switch len(tabs) {
+ case 0:
+ return new(Table)
+
+ case 1:
+ return tabs[0]
+ }
+
+ // Construct each column.
+ var out Builder
+ seqs := make([]slice.T, len(tabs))
+ for _, col := range tabs[0].Columns() {
+ seqs = seqs[:0]
+ for _, tab := range tabs {
+ seqs = append(seqs, tab.Column(col))
+ }
+ out.Add(col, slice.Concat(seqs...))
+ }
+
+ return out.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/head.go b/vendor/github.com/aclements/go-gg/table/head.go
new file mode 100644
index 0000000..3c1b2ec
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/head.go
@@ -0,0 +1,69 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import "reflect"
+
+// Head returns the first n rows in each Table of g.
+func Head(g Grouping, n int) Grouping {
+ return headTail(g, n, false)
+}
+
+// Tail returns the last n rows in each Table of g.
+func Tail(g Grouping, n int) Grouping {
+ return headTail(g, n, true)
+}
+
+func headTail(g Grouping, n int, tail bool) Grouping {
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ if t.Len() <= n {
+ return t
+ }
+
+ var nt Builder
+ for _, col := range t.Columns() {
+ if cv, ok := t.Const(col); ok {
+ nt.AddConst(col, cv)
+ continue
+ }
+
+ cv := reflect.ValueOf(t.Column(col))
+ if tail {
+ cv = cv.Slice(t.Len()-n, t.Len())
+ } else {
+ cv = cv.Slice(0, n)
+ }
+ nt.Add(col, cv.Interface())
+ }
+ return nt.Done()
+ })
+}
+
+// HeadTables returns the first n tables in g.
+func HeadTables(g Grouping, n int) Grouping {
+ return headTailTables(g, n, false)
+}
+
+// TailTables returns the first n tables in g.
+func TailTables(g Grouping, n int) Grouping {
+ return headTailTables(g, n, true)
+}
+
+func headTailTables(g Grouping, n int, tail bool) Grouping {
+ tables := g.Tables()
+ if len(tables) <= n {
+ return g
+ } else if tail {
+ tables = tables[len(tables)-n:]
+ } else {
+ tables = tables[:n]
+ }
+
+ var ng GroupingBuilder
+ for _, gid := range tables {
+ ng.Add(gid, g.Table(gid))
+ }
+ return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/join.go b/vendor/github.com/aclements/go-gg/table/join.go
new file mode 100644
index 0000000..03ce51b
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/join.go
@@ -0,0 +1,77 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+// Join joins g1 and g2 on tables with identical group IDs where col1
+// in g1 equals col2 in g2. It maintains the group order of g1, except
+// that groups that aren't in g2 are removed, and maintains the row
+// order of g1, followed by the row order of g2.
+//
+// TODO: Support join on more than one column.
+func Join(g1 Grouping, col1 string, g2 Grouping, col2 string) Grouping {
+ var ng GroupingBuilder
+ for _, gid := range g1.Tables() {
+ t1, t2 := g1.Table(gid), g2.Table(gid)
+ if t2 == nil {
+ continue
+ }
+
+ // TODO: Optimize for cases where col1 and/or col2 are
+ // constant.
+
+ // Index col2 in t2.
+ ridx := make(map[interface{}][]int)
+ rv := reflect.ValueOf(t2.MustColumn(col2))
+ for i, l := 0, rv.Len(); i < l; i++ {
+ v := rv.Index(i).Interface()
+ ridx[v] = append(ridx[v], i)
+ }
+
+ // For each row in t1, find the matching rows in col2
+ // and build up the row indexes for t1 and t2.
+ idx1, idx2 := []int{}, []int{}
+ lv := reflect.ValueOf(t1.MustColumn(col1))
+ for i, l := 0, lv.Len(); i < l; i++ {
+ r := ridx[lv.Index(i).Interface()]
+ for range r {
+ idx1 = append(idx1, i)
+ }
+ idx2 = append(idx2, r...)
+ }
+
+ // Build the joined table.
+ var nt Builder
+ for _, col := range t1.Columns() {
+ if cv, ok := t1.Const(col); ok {
+ nt.Add(col, cv)
+ continue
+ }
+ nt.Add(col, slice.Select(t1.Column(col), idx1))
+ }
+ for _, col := range t2.Columns() {
+ // Often the join column is the same in both
+ // and we can skip it because we added it from
+ // the first table.
+ if col == col1 && col == col2 {
+ continue
+ }
+
+ if cv, ok := t2.Const(col); ok {
+ nt.Add(col, cv)
+ continue
+ }
+ nt.Add(col, slice.Select(t2.Column(col), idx2))
+ }
+
+ ng.Add(gid, nt.Done())
+ }
+ return ng.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/map.go b/vendor/github.com/aclements/go-gg/table/map.go
new file mode 100644
index 0000000..87641a5
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/map.go
@@ -0,0 +1,170 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "fmt"
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// MapTables applies f to each Table in g and returns a new Grouping
+// with the same group structure as g, but with the Tables returned by
+// f.
+func MapTables(g Grouping, f func(gid GroupID, table *Table) *Table) Grouping {
+ var out GroupingBuilder
+ for _, gid := range g.Tables() {
+ out.Add(gid, f(gid, g.Table(gid)))
+ }
+ return out.Done()
+}
+
+// MapCols applies f to a set of input columns to construct a set of
+// new output columns.
+//
+// For each Table in g, MapCols calls f(in[0], in[1], ..., out[0],
+// out[1], ...) where in[i] is column incols[i]. f should process the
+// values in the input column slices and fill output columns slices
+// out[j] accordingly. MapCols returns a new Grouping that adds each
+// outcols[j] bound to out[j].
+//
+// If all of the input columns are constant for a given table, MapCols
+// will call f with all slices of length 1. The input column slices
+// will contain the constant column values and MapCols will bind each
+// output column value out[i][0] as a constant.
+func MapCols(g Grouping, f interface{}, incols ...string) func(outcols ...string) Grouping {
+ return func(outcols ...string) Grouping {
+ fv := reflect.ValueOf(f)
+ if fv.Kind() != reflect.Func {
+ panic(&generic.TypeError{fv.Type(), nil, "must be a function"})
+ }
+ ft := fv.Type()
+ if ft.NumIn() != len(incols)+len(outcols) {
+ panic(&generic.TypeError{ft, nil, fmt.Sprintf("has the wrong number of arguments; expected %d", len(incols)+len(outcols))})
+ }
+ if ft.NumOut() != 0 {
+ panic(&generic.TypeError{ft, nil, "has the wrong number of results; expected 0"})
+ }
+
+ // Create output column slices.
+ totalRows := 0
+ for _, gid := range g.Tables() {
+ t := g.Table(gid)
+ colloop:
+ for _, incol := range incols {
+ if _, ok := t.Const(incol); !ok {
+ totalRows += g.Table(gid).Len()
+ break colloop
+ }
+ }
+ }
+ ocols := make([]reflect.Value, len(outcols))
+ for i := range ocols {
+ ocols[i] = reflect.MakeSlice(ft.In(i+len(incols)), totalRows, totalRows)
+ }
+
+ // Apply f to each group.
+ var out GroupingBuilder
+ args := make([]reflect.Value, len(incols)+len(outcols))
+ opos := 0
+ for _, gid := range g.Tables() {
+ t := g.Table(gid)
+
+ // Are all inputs are constants?
+ allConst := true
+ for _, incol := range incols {
+ if _, ok := t.Const(incol); !ok {
+ allConst = false
+ break
+ }
+ }
+ if allConst {
+ for i, incol := range incols {
+ cv, _ := t.Const(incol)
+ args[i] = reflect.MakeSlice(ColType(t, incol), 1, 1)
+ args[i].Index(0).Set(reflect.ValueOf(cv))
+ }
+ for i, ocol := range ocols {
+ args[i+len(incols)] = reflect.MakeSlice(ocol.Type(), 1, 1)
+ }
+
+ fv.Call(args)
+
+ tb := NewBuilder(t)
+ for i, outcol := range outcols {
+ tb.AddConst(outcol, args[i+len(incols)].Index(0).Interface())
+ }
+ out.Add(gid, tb.Done())
+ continue
+ }
+
+ // Prepare arguments.
+ for i, incol := range incols {
+ args[i] = reflect.ValueOf(t.MustColumn(incol))
+ }
+ for i, ocol := range ocols {
+ args[i+len(incols)] = ocol.Slice(opos, opos+t.Len())
+ }
+ opos += t.Len()
+
+ // Call f.
+ fv.Call(args)
+
+ // Add output columns.
+ tb := NewBuilder(t)
+ for i, outcol := range outcols {
+ tb.Add(outcol, args[i+len(incols)].Interface())
+ }
+ out.Add(gid, tb.Done())
+ }
+ return out.Done()
+ }
+}
+
+// Rename returns g with column 'from' renamed to 'to'. The column
+// retains its position.
+func Rename(g Grouping, from, to string) Grouping {
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ t.MustColumn(from)
+ var nt Builder
+ for _, col := range t.Columns() {
+ if col == to {
+ continue
+ }
+
+ ncol := col
+ if col == from {
+ ncol = to
+ }
+
+ if cv, ok := t.Const(col); ok {
+ nt.AddConst(ncol, cv)
+ } else {
+ nt.Add(ncol, t.Column(col))
+ }
+ }
+ return nt.Done()
+ })
+}
+
+// Remove returns g with column 'col' removed.
+func Remove(g Grouping, col string) Grouping {
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ t.MustColumn(col)
+ var nt Builder
+ for _, col2 := range t.Columns() {
+ if col == col2 {
+ continue
+ }
+ if cv, ok := t.Const(col2); ok {
+ nt.AddConst(col2, cv)
+ } else {
+ nt.Add(col2, t.Column(col2))
+ }
+ }
+ return nt.Done()
+ })
+}
diff --git a/vendor/github.com/aclements/go-gg/table/new.go b/vendor/github.com/aclements/go-gg/table/new.go
new file mode 100644
index 0000000..8094593
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/new.go
@@ -0,0 +1,105 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "reflect"
+ "strconv"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// TableFromStructs converts a []T where T is a struct to a Table
+// where the columns of the table correspond to T's exported fields.
+func TableFromStructs(structs Slice) *Table {
+ s := reflectSlice(structs)
+ st := s.Type()
+ if st.Elem().Kind() != reflect.Struct {
+ panic(&generic.TypeError{st, nil, "is not a slice of struct"})
+ }
+
+ var t Builder
+ rows := s.Len()
+ var rec func(reflect.Type, []int)
+ rec = func(typ reflect.Type, index []int) {
+ for fn := 0; fn < typ.NumField(); fn++ {
+ field := typ.Field(fn)
+ if field.PkgPath != "" {
+ continue
+ }
+ oldIndexLen := len(index)
+ index = append(index, field.Index...)
+ if field.Anonymous {
+ rec(field.Type, index)
+ } else {
+ col := reflect.MakeSlice(reflect.SliceOf(field.Type), rows, rows)
+ for i := 0; i < rows; i++ {
+ col.Index(i).Set(s.Index(i).FieldByIndex(index))
+ }
+ t.Add(field.Name, col.Interface())
+ }
+ index = index[:oldIndexLen]
+ }
+ }
+ rec(st.Elem(), []int{})
+ return t.Done()
+}
+
+// TableFromStrings converts a [][]string to a Table. This is intended
+// for processing external data, such as from CSV files. If coerce is
+// true, TableFromStrings will convert columns to []int or []float
+// when every string in that column is accepted by strconv.ParseInt or
+// strconv.ParseFloat, respectively.
+func TableFromStrings(cols []string, rows [][]string, coerce bool) *Table {
+ var t Builder
+ for i, col := range cols {
+ slice := make([]string, len(rows))
+ for j, row := range rows {
+ slice[j] = row[i]
+ }
+
+ var colData interface{} = slice
+ switch {
+ case coerce && len(slice) > 0:
+ // Try []int.
+ var err error
+ for _, str := range slice {
+ _, err = strconv.ParseInt(str, 10, 0)
+ if err != nil {
+ break
+ }
+ }
+ if err == nil {
+ nslice := make([]int, len(rows))
+ for i, str := range slice {
+ v, _ := strconv.ParseInt(str, 10, 0)
+ nslice[i] = int(v)
+ }
+ colData = nslice
+ break
+ }
+
+ // Try []float64. This must be done after
+ // []int. It's also more expensive.
+ for _, str := range slice {
+ _, err = strconv.ParseFloat(str, 64)
+ if err != nil {
+ break
+ }
+ }
+ if err == nil {
+ nslice := make([]float64, len(rows))
+ for i, str := range slice {
+ nslice[i], _ = strconv.ParseFloat(str, 64)
+ }
+ colData = nslice
+ break
+ }
+ }
+
+ t.Add(col, colData)
+ }
+ return t.Done()
+}
diff --git a/vendor/github.com/aclements/go-gg/table/pivot.go b/vendor/github.com/aclements/go-gg/table/pivot.go
new file mode 100644
index 0000000..6f488f9
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/pivot.go
@@ -0,0 +1,153 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+)
+
+// Pivot converts rows of g into columns. label and value must name
+// columns in g, and the label column must have type []string. Pivot
+// returns a Grouping with a new column named after each distinct
+// value in the label column, where the values in that column
+// correspond to the values from the value column. All other columns
+// (besides label and value) are copied to the output. If, for a given
+// column in an output row, no input row has that column in the label
+// column, the output cell will have the zero value for its type.
+func Pivot(g Grouping, label, value string) Grouping {
+ // Find all unique values of label. These are the new columns.
+ labels := []string{}
+ lset := map[string]int{}
+ for _, gid := range g.Tables() {
+ for _, l := range g.Table(gid).MustColumn(label).([]string) {
+ if _, ok := lset[l]; !ok {
+ lset[l] = len(lset)
+ labels = append(labels, l)
+ }
+ }
+ }
+
+ // Get all columns that are not label or value.
+ groupCols := []string{}
+ for _, col := range g.Columns() {
+ if col != label && col != value {
+ groupCols = append(groupCols, col)
+ }
+ }
+
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ var nt Builder
+
+ // Group by all other columns. Each group in gg
+ // becomes an output row.
+ gg := GroupBy(t, groupCols...)
+
+ // Copy grouped-by values.
+ for _, groupCol := range groupCols {
+ cv := reflect.MakeSlice(reflect.TypeOf(t.Column(groupCol)), len(gg.Tables()), len(gg.Tables()))
+ for i, gid := range gg.Tables() {
+ sub := gg.Table(gid)
+ cv.Index(i).Set(reflect.ValueOf(sub.Column(groupCol)).Index(0))
+ }
+ nt.Add(groupCol, cv.Interface())
+ }
+
+ // Initialize new columns.
+ newCols := make([]reflect.Value, len(lset))
+ vt := reflect.TypeOf(t.MustColumn(value))
+ for i := range newCols {
+ newCols[i] = reflect.MakeSlice(vt, len(gg.Tables()), len(gg.Tables()))
+ }
+
+ // Fill in new columns.
+ for i, gid := range gg.Tables() {
+ sub := gg.Table(gid)
+
+ vcol := reflect.ValueOf(sub.MustColumn(value))
+ for j, l := range sub.MustColumn(label).([]string) {
+ val := vcol.Index(j)
+ newCols[lset[l]].Index(i).Set(val)
+ }
+ }
+
+ // Add new columns to output table.
+ for i, newCol := range newCols {
+ nt.Add(labels[i], newCol.Interface())
+ }
+
+ return nt.Done()
+ })
+}
+
+// Unpivot converts columns of g into rows. The returned Grouping
+// consists of the columns of g *not* listed in cols, plus two columns
+// named by the label and value arguments. For each input row in g,
+// the returned Grouping will have len(cols) output rows. The i'th
+// such output row corresponds to column cols[i] in the input row. The
+// label column will contain the name of the unpivoted column,
+// cols[i], and the value column will contain that column's value from
+// the input row. The values of all other columns in the input row
+// will be repeated across the output rows. All columns in cols must
+// have the same type.
+func Unpivot(g Grouping, label, value string, cols ...string) Grouping {
+ if len(cols) == 0 {
+ panic("Unpivot requires at least 1 column")
+ }
+
+ colSet := map[string]bool{}
+ for _, col := range cols {
+ colSet[col] = true
+ }
+
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ var nt Builder
+
+ // Repeat all other columns len(cols) times.
+ ntlen := t.Len() * len(cols)
+ for _, name := range t.Columns() {
+ if colSet[name] || name == label || name == value {
+ continue
+ }
+
+ col := reflect.ValueOf(t.Column(name))
+ ncol := reflect.MakeSlice(col.Type(), ntlen, ntlen)
+ for i, l := 0, col.Len(); i < l; i++ {
+ v := col.Index(i)
+ for j := range cols {
+ ncol.Index(i*len(cols) + j).Set(v)
+ }
+ }
+
+ nt.Add(name, ncol.Interface())
+ }
+
+ // Get input columns.
+ var vt reflect.Type
+ colvs := make([]reflect.Value, len(cols))
+ for i, col := range cols {
+ colvs[i] = reflect.ValueOf(t.MustColumn(col))
+ if i == 0 {
+ vt = colvs[i].Type()
+ } else if vt != colvs[i].Type() {
+ panic(&generic.TypeError{vt, colvs[i].Type(), "; cannot Unpivot columns with different types"})
+ }
+ }
+
+ // Create label and value columns.
+ lcol := make([]string, 0, ntlen)
+ vcol := reflect.MakeSlice(vt, ntlen, ntlen)
+ for i := 0; i < t.Len(); i++ {
+ lcol = append(lcol, cols...)
+ for j, colv := range colvs {
+ vcol.Index(i*len(cols) + j).Set(colv.Index(i))
+ }
+ }
+ nt.Add(label, lcol).Add(value, vcol.Interface())
+
+ return nt.Done()
+ })
+}
diff --git a/vendor/github.com/aclements/go-gg/table/print.go b/vendor/github.com/aclements/go-gg/table/print.go
new file mode 100644
index 0000000..1795865
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/print.go
@@ -0,0 +1,119 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "fmt"
+ "io"
+ "os"
+ "reflect"
+ "strings"
+)
+
+// TODO: Have a format struct with options for things like column
+// separator, and header separator. Provide some defaults ones for,
+// e.g., Markdown, CSV, TSV, and such. Make top-level Print and Fprint
+// call methods in some default format.
+
+// Print(...) is shorthand for Fprint(os.Stderr, ...).
+func Print(g Grouping, formats ...string) error {
+ return Fprint(os.Stdout, g, formats...)
+}
+
+// Fprint prints Grouping g to w. formats[i] specifies a fmt-style
+// format string for column i. If there are more columns than formats,
+// remaining columns are formatted with %v (in particular, formats may
+// be omitted entirely to use %v for all columns). Numeric columns are
+// right aligned; all other column types are left aligned.
+func Fprint(w io.Writer, g Grouping, formats ...string) error {
+ if g.Columns() == nil {
+ return nil
+ }
+
+ // Convert each column to strings.
+ ss := make([][]string, len(g.Columns()))
+ rowFmts := make([]string, len(g.Columns()))
+ for i, col := range g.Columns() {
+ format := "%v"
+ if i < len(formats) {
+ format = formats[i]
+ }
+
+ // Format column.
+ var valKind reflect.Kind
+ ss[i] = []string{col}
+ for _, gid := range g.Tables() {
+ seq := reflect.ValueOf(g.Table(gid).Column(col))
+ for row := 0; row < seq.Len(); row++ {
+ str := fmt.Sprintf(format, seq.Index(row).Interface())
+ ss[i] = append(ss[i], str)
+ }
+
+ if valKind == reflect.Invalid {
+ valKind = seq.Type().Elem().Kind()
+ }
+ }
+
+ // Find column width.
+ width := 0
+ for _, s := range ss[i] {
+ if len(s) > width {
+ width = len(s)
+ }
+ }
+
+ // If it's a numeric column, right align.
+ //
+ // TODO: Even better would be to decimal align, though
+ // that may require some understanding of the format;
+ // or we could only do it for the default format.
+ switch valKind {
+ case reflect.Int, reflect.Int8, reflect.Int16, reflect.Int32, reflect.Int64, reflect.Uint, reflect.Uint8, reflect.Uint16, reflect.Uint32, reflect.Uint64, reflect.Uintptr, reflect.Float32, reflect.Float64:
+ width = -width
+ }
+
+ if i == len(g.Columns())-1 && width > 0 {
+ // Don't pad the last column.
+ rowFmts[i] = "%s"
+ } else {
+ rowFmts[i] = fmt.Sprintf("%%%ds", -width)
+ }
+ }
+
+ // Compute group headers.
+ groups := []GroupID{}
+ groupPos := []int{}
+ lastPos := 1
+ for _, gid := range g.Tables() {
+ groups = append(groups, gid)
+ groupPos = append(groupPos, lastPos)
+ lastPos += g.Table(gid).Len()
+ }
+ if len(groups) == 1 && groups[0] == RootGroupID {
+ groups, groupPos = nil, nil
+ }
+
+ // Print rows.
+ rowFmt := strings.Join(rowFmts, " ") + "\n"
+ rowBuf := make([]interface{}, len(rowFmts))
+ for row := 0; row < len(ss[0]); row++ {
+ if len(groupPos) > 0 && row == groupPos[0] {
+ _, err := fmt.Fprintf(w, "-- %s\n", groups[0])
+ if err != nil {
+ return err
+ }
+ groups, groupPos = groups[1:], groupPos[1:]
+ }
+
+ for col := range rowBuf {
+ rowBuf[col] = ss[col][row]
+ }
+ _, err := fmt.Fprintf(w, rowFmt, rowBuf...)
+ if err != nil {
+ return err
+ }
+ }
+ return nil
+}
diff --git a/vendor/github.com/aclements/go-gg/table/sort.go b/vendor/github.com/aclements/go-gg/table/sort.go
new file mode 100644
index 0000000..ec64e00
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/sort.go
@@ -0,0 +1,93 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package table
+
+import (
+ "sort"
+
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+// SortBy sorts each group of g by the named columns. If a column's
+// type implements sort.Interface, rows will be sorted according to
+// that order. Otherwise, the values in the column must be naturally
+// ordered (their types must be orderable by the Go specification). If
+// neither is true, SortBy panics with a *generic.TypeError. If more
+// than one column is given, SortBy sorts by the tuple of the columns;
+// that is, if two values in the first column are equal, they are
+// sorted by the second column, and so on.
+func SortBy(g Grouping, cols ...string) Grouping {
+ // Sort each group.
+ sorters := make([]sort.Interface, len(cols))
+ return MapTables(g, func(_ GroupID, t *Table) *Table {
+ // Create sorters for each column.
+ sorters = sorters[:0]
+ for _, col := range cols {
+ if _, ok := t.Const(col); ok {
+ continue
+ }
+ seq := t.MustColumn(col)
+ sorter := slice.Sorter(seq)
+ if sort.IsSorted(sorter) {
+ continue
+ }
+ sorters = append(sorters, sorter)
+ }
+
+ if len(sorters) == 0 {
+ // Avoid shuffling everything by the identity
+ // permutation.
+ return t
+ }
+
+ // Generate an initial permutation sequence.
+ perm := make([]int, t.Len())
+ for i := range perm {
+ perm[i] = i
+ }
+
+ // Sort the permutation sequence.
+ sort.Stable(&permSort{perm, sorters})
+
+ // Permute all columns.
+ var nt Builder
+ for _, name := range t.Columns() {
+ if cv, ok := t.Const(name); ok {
+ nt.AddConst(name, cv)
+ continue
+ }
+ seq := t.Column(name)
+ seq = slice.Select(seq, perm)
+ nt.Add(name, seq)
+ }
+ return nt.Done()
+ })
+}
+
+type permSort struct {
+ perm []int
+ keys []sort.Interface
+}
+
+func (s *permSort) Len() int {
+ return len(s.perm)
+}
+
+func (s *permSort) Less(i, j int) bool {
+ // Since there's no way to ask about equality, we have to do
+ // extra work for all of the keys except the last.
+ for _, key := range s.keys[:len(s.keys)-1] {
+ if key.Less(s.perm[i], s.perm[j]) {
+ return true
+ } else if key.Less(s.perm[j], s.perm[i]) {
+ return false
+ }
+ }
+ return s.keys[len(s.keys)-1].Less(s.perm[i], s.perm[j])
+}
+
+func (s *permSort) Swap(i, j int) {
+ s.perm[i], s.perm[j] = s.perm[j], s.perm[i]
+}
diff --git a/vendor/github.com/aclements/go-gg/table/table.go b/vendor/github.com/aclements/go-gg/table/table.go
new file mode 100644
index 0000000..5821752
--- /dev/null
+++ b/vendor/github.com/aclements/go-gg/table/table.go
@@ -0,0 +1,489 @@
+// Copyright 2016 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package table implements ordered, grouped two dimensional relations.
+//
+// There are two related abstractions: Table and Grouping.
+//
+// A Table is an ordered relation of rows and columns. Each column is
+// a Go slice and hence must be homogeneously typed, but different
+// columns may have different types. All columns in a Table have the
+// same number of rows.
+//
+// A Grouping generalizes a Table by grouping the Table's rows into
+// zero or more groups. A Table is itself a Grouping with zero or one
+// groups. Most operations take a Grouping and operate on each group
+// independently, though some operations sub-divide or combine groups.
+//
+// The structures of both Tables and Groupings are immutable. They are
+// constructed using a Builder or a GroupingBuilder, respectively, and
+// then "frozen" into their respective immutable data structures.
+package table
+
+import (
+ "fmt"
+ "reflect"
+
+ "github.com/aclements/go-gg/generic"
+ "github.com/aclements/go-gg/generic/slice"
+)
+
+// TODO
+//
+// Rename Table to T?
+//
+// Make Table an interface? Then columns could be constructed lazily.
+//
+// Do all transformation functions as func(g Grouping) Grouping? That
+// could be a "Transform" type that has easy methods for chaining. In
+// a lot of cases, transformation functions could just return the
+// Transform returned by another function (like MapTables).
+//
+// Make an error type for "unknown column".
+
+// A Table is an immutable, ordered two dimensional relation. It
+// consists of a set of named columns. Each column is a sequence of
+// values of a consistent type or a constant value. All (non-constant)
+// columns have the same length.
+//
+// The zero value of Table is the "empty table": it has no rows and no
+// columns. Note that a Table may have one or more columns, but no
+// rows; such a Table is *not* considered empty.
+//
+// A Table is also a trivial Grouping. If a Table is empty, it has no
+// groups and hence the zero value of Table is also the "empty group".
+// Otherwise, it consists only of the root group, RootGroupID.
+type Table struct {
+ cols map[string]Slice
+ consts map[string]interface{}
+ colNames []string
+ len int
+}
+
+// A Builder constructs a Table one column at a time.
+//
+// The zero value of a Builder represents an empty Table.
+type Builder struct {
+ t Table
+}
+
+// A Grouping is an immutable set of tables with identical sets of
+// columns, each identified by a distinct GroupID.
+//
+// Visually, a Grouping can be thought of as follows:
+//
+// Col A Col B Col C
+// ------ group /a ------
+// 0 5.4 "x" 90
+// 1 -.2 "y" 30
+// ------ group /b ------
+// 0 9.3 "a" 10
+//
+// Like a Table, a Grouping's structure is immutable. To construct a
+// Grouping, use a GroupingBuilder.
+//
+// Despite the fact that GroupIDs form a hierarchy, a Grouping ignores
+// this hierarchy and simply operates on a flat map of distinct
+// GroupIDs to Tables.
+type Grouping interface {
+ // Columns returns the names of the columns in this Grouping,
+ // or nil if there are no Tables or the group consists solely
+ // of empty Tables. All Tables in this Grouping have the same
+ // set of columns.
+ Columns() []string
+
+ // Tables returns the group IDs of the tables in this
+ // Grouping.
+ Tables() []GroupID
+
+ // Table returns the Table in group gid, or nil if there is no
+ // such Table.
+ Table(gid GroupID) *Table
+}
+
+// A GroupingBuilder constructs a Grouping one table a time.
+//
+// The zero value of a GroupingBuilder represents an empty Grouping
+// with no tables and no columns.
+type GroupingBuilder struct {
+ g groupedTable
+ colTypes []reflect.Type
+}
+
+type groupedTable struct {
+ tables map[GroupID]*Table
+ groups []GroupID
+ colNames []string
+}
+
+// A Slice is a Go slice value.
+//
+// This is primarily for documentation. There is no way to statically
+// enforce this in Go; however, functions that expect a Slice will
+// panic with a *generic.TypeError if passed a non-slice value.
+type Slice interface{}
+
+func reflectSlice(s Slice) reflect.Value {
+ rv := reflect.ValueOf(s)
+ if rv.Kind() != reflect.Slice {
+ panic(&generic.TypeError{rv.Type(), nil, "is not a slice"})
+ }
+ return rv
+}
+
+// NewBuilder returns a new Builder. If t is non-nil, it populates the
+// new Builder with the columns from t.
+func NewBuilder(t *Table) *Builder {
+ if t == nil {
+ return new(Builder)
+ }
+ b := Builder{Table{
+ cols: make(map[string]Slice),
+ consts: make(map[string]interface{}),
+ colNames: append([]string(nil), t.Columns()...),
+ len: t.len,
+ }}
+ for k, v := range t.cols {
+ b.t.cols[k] = v
+ }
+ for k, v := range t.consts {
+ b.t.consts[k] = v
+ }
+ return &b
+}
+
+// Add adds a column to b, or removes the named column if data is nil.
+// If b already has a column with the given name, Add replaces it. If
+// data is non-nil, it must have the same length as any existing
+// columns or Add will panic.
+func (b *Builder) Add(name string, data Slice) *Builder {
+ if data == nil {
+ // Remove the column.
+ if _, ok := b.t.cols[name]; !ok {
+ if _, ok := b.t.consts[name]; !ok {
+ // Nothing to remove.
+ return b
+ }
+ }
+ delete(b.t.cols, name)
+ delete(b.t.consts, name)
+ for i, n := range b.t.colNames {
+ if n == name {
+ copy(b.t.colNames[i:], b.t.colNames[i+1:])
+ b.t.colNames = b.t.colNames[:len(b.t.colNames)-1]
+ break
+ }
+ }
+ return b
+ }
+
+ // Are we replacing an existing column?
+ _, replace := b.t.cols[name]
+ if !replace {
+ _, replace = b.t.consts[name]
+ }
+
+ // Check the column and add it.
+ rv := reflectSlice(data)
+ dataLen := rv.Len()
+ if len(b.t.cols) == 0 || (replace && len(b.t.cols) == 1) {
+ if b.t.cols == nil {
+ b.t.cols = make(map[string]Slice)
+ }
+ // First non-constant column (possibly replacing the
+ // only non-constant column).
+ b.t.cols[name] = data
+ b.t.len = dataLen
+ } else if b.t.len != dataLen {
+ panic(fmt.Sprintf("cannot add column %q with %d elements to table with %d rows", name, dataLen, b.t.len))
+ } else {
+ b.t.cols[name] = data
+ }
+
+ if replace {
+ // Make sure it's not in constants.
+ delete(b.t.consts, name)
+ } else {
+ b.t.colNames = append(b.t.colNames, name)
+ }
+
+ return b
+}
+
+// AddConst adds a constant column to b whose value is val. If b
+// already has a column with this name, AddConst replaces it.
+//
+// A constant column has the same value in every row of the Table. It
+// does not itself have an inherent length.
+func (b *Builder) AddConst(name string, val interface{}) *Builder {
+ // Are we replacing an existing column?
+ _, replace := b.t.cols[name]
+ if !replace {
+ _, replace = b.t.consts[name]
+ }
+
+ if b.t.consts == nil {
+ b.t.consts = make(map[string]interface{})
+ }
+ b.t.consts[name] = val
+
+ if replace {
+ // Make sure it's not in cols.
+ delete(b.t.cols, name)
+ } else {
+ b.t.colNames = append(b.t.colNames, name)
+ }
+
+ return b
+}
+
+// Has returns true if b has a column named "name".
+func (b *Builder) Has(name string) bool {
+ if _, ok := b.t.cols[name]; ok {
+ return true
+ }
+ if _, ok := b.t.consts[name]; ok {
+ return true
+ }
+ return false
+}
+
+// Done returns the constructed Table and resets b.
+func (b *Builder) Done() *Table {
+ if len(b.t.colNames) == 0 {
+ return new(Table)
+ }
+ t := b.t
+ b.t = Table{}
+ return &t
+}
+
+// Len returns the number of rows in Table t.
+func (t *Table) Len() int {
+ return t.len
+}
+
+// Columns returns the names of the columns in Table t, or nil if this
+// Table is empty.
+func (t *Table) Columns() []string {
+ return t.colNames
+}
+
+// Column returns the slice of data in column name of Table t, or nil
+// if there is no such column. If name is a constant column, this
+// returns a slice with the constant value repeated to the length of
+// the Table.
+func (t *Table) Column(name string) Slice {
+ if c, ok := t.cols[name]; ok {
+ // It's a regular column or a constant column with a
+ // cached expansion.
+ return c
+ }
+
+ if cv, ok := t.consts[name]; ok {
+ // Expand the constant column and cache the result.
+ expanded := slice.Repeat(cv, t.len)
+ t.cols[name] = expanded
+ return expanded
+ }
+
+ return nil
+}
+
+// MustColumn is like Column, but panics if there is no such column.
+func (t *Table) MustColumn(name string) Slice {
+ if c := t.Column(name); c != nil {
+ return c
+ }
+ panic(fmt.Sprintf("unknown column %q", name))
+}
+
+// Const returns the value of constant column name. If this column
+// does not exist or is not a constant column, Const returns nil,
+// false.
+func (t *Table) Const(name string) (val interface{}, ok bool) {
+ cv, ok := t.consts[name]
+ return cv, ok
+}
+
+// isEmpty returns true if t is an empty Table, meaning it has no rows
+// or columns.
+func (t *Table) isEmpty() bool {
+ return t.colNames == nil
+}
+
+// Tables returns the groups IDs in this Table. If t is empty, there
+// are no group IDs. Otherwise, there is only RootGroupID.
+func (t *Table) Tables() []GroupID {
+ if t.isEmpty() {
+ return []GroupID{}
+ }
+ return []GroupID{RootGroupID}
+}
+
+// Table returns t if gid is RootGroupID and t is not empty; otherwise
+// it returns nil.
+func (t *Table) Table(gid GroupID) *Table {
+ if gid == RootGroupID && !t.isEmpty() {
+ return t
+ }
+ return nil
+}
+
+// NewGroupingBuilder returns a new GroupingBuilder. If g is non-nil,
+// it populates the new GroupingBuilder with the tables from g.
+func NewGroupingBuilder(g Grouping) *GroupingBuilder {
+ if g == nil {
+ return new(GroupingBuilder)
+ }
+ b := GroupingBuilder{groupedTable{
+ tables: make(map[GroupID]*Table),
+ groups: append([]GroupID(nil), g.Tables()...),
+ colNames: append([]string(nil), g.Columns()...),
+ }, nil}
+ for _, gid := range g.Tables() {
+ t := g.Table(gid)
+ b.g.tables[gid] = t
+ if b.colTypes == nil {
+ b.colTypes = colTypes(t)
+ }
+ }
+ return &b
+}
+
+func colTypes(t *Table) []reflect.Type {
+ colTypes := make([]reflect.Type, len(t.colNames))
+ for i, col := range t.colNames {
+ if c, ok := t.cols[col]; ok {
+ colTypes[i] = reflect.TypeOf(c).Elem()
+ } else {
+ colTypes[i] = reflect.TypeOf(t.consts[col])
+ }
+ }
+ return colTypes
+}
+
+// Add adds a Table to b, or removes a table if t is nil. If t is the
+// empty Table, this is a no-op because the empty Table contains no
+// groups. If gid already exists, Add replaces it. Table t must have
+// the same columns as any existing Tables in this Grouping and they
+// must have identical types; otherwise, Add will panic.
+//
+// TODO This doesn't make it easy to combine two Groupings. It could
+// instead take a Grouping and reparent it.
+func (b *GroupingBuilder) Add(gid GroupID, t *Table) *GroupingBuilder {
+ if t == nil {
+ if _, ok := b.g.tables[gid]; !ok {
+ // Nothing to remove.
+ return b
+ }
+ delete(b.g.tables, gid)
+ for i, g2 := range b.g.groups {
+ if g2 == gid {
+ copy(b.g.groups[i:], b.g.groups[i+1:])
+ b.g.groups = b.g.groups[:len(b.g.groups)-1]
+ break
+ }
+ }
+ return b
+ }
+
+ if t != nil && t.isEmpty() {
+ // Adding an empty table has no effect.
+ return b
+ }
+
+ if len(b.g.groups) == 1 && b.g.groups[0] == gid {
+ // We're replacing the only group. This is allowed to
+ // change the shape of the Grouping.
+ b.g.tables[gid] = t
+ b.g.colNames = t.Columns()
+ b.colTypes = colTypes(t)
+ return b
+ } else if len(b.g.groups) == 0 {
+ b.g.tables = map[GroupID]*Table{gid: t}
+ b.g.groups = []GroupID{gid}
+ b.g.colNames = t.Columns()
+ b.colTypes = colTypes(t)
+ return b
+ }
+
+ // Check that t's column names match.
+ matches := true
+ if len(t.colNames) != len(b.g.colNames) {
+ matches = false
+ } else {
+ for i, n := range t.colNames {
+ if b.g.colNames[i] != n {
+ matches = false
+ break
+ }
+ }
+ }
+ if !matches {
+ panic(fmt.Sprintf("table columns %q do not match group columns %q", t.colNames, b.g.colNames))
+ }
+
+ // Check that t's column types match.
+ for i, col := range b.g.colNames {
+ t0 := b.colTypes[i]
+ var t1 reflect.Type
+ if c, ok := t.cols[col]; ok {
+ t1 = reflect.TypeOf(c).Elem()
+ } else if cv, ok := t.consts[col]; ok {
+ t1 = reflect.TypeOf(cv)
+ }
+ if t0 != t1 {
+ panic(&generic.TypeError{t0, t1, fmt.Sprintf("for column %q are not the same", col)})
+ }
+ }
+
+ // Add t.
+ if _, ok := b.g.tables[gid]; !ok {
+ b.g.groups = append(b.g.groups, gid)
+ }
+ b.g.tables[gid] = t
+
+ return b
+}
+
+// Done returns the constructed Grouping and resets b.
+func (b *GroupingBuilder) Done() Grouping {
+ if len(b.g.groups) == 0 {
+ return new(groupedTable)
+ }
+ g := b.g
+ b.g = groupedTable{}
+ return &g
+}
+
+func (g *groupedTable) Columns() []string {
+ return g.colNames
+}
+
+func (g *groupedTable) Tables() []GroupID {
+ return g.groups
+}
+
+func (g *groupedTable) Table(gid GroupID) *Table {
+ return g.tables[gid]
+}
+
+// ColType returns the type of column col in g. This will always be a
+// slice type, even if col is a constant column. ColType panics if col
+// is unknown.
+//
+// TODO: If I introduce a first-class representation for a grouped
+// column, this should probably be in that.
+func ColType(g Grouping, col string) reflect.Type {
+ tables := g.Tables()
+ if len(tables) == 0 {
+ panic(fmt.Sprintf("unknown column %q", col))
+ }
+ t0 := g.Table(tables[0])
+ if cv, ok := t0.Const(col); ok {
+ return reflect.SliceOf(reflect.TypeOf(cv))
+ }
+ return reflect.TypeOf(t0.MustColumn(col))
+}
diff --git a/vendor/github.com/aclements/go-moremath/LICENSE b/vendor/github.com/aclements/go-moremath/LICENSE
new file mode 100644
index 0000000..d29b372
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/LICENSE
@@ -0,0 +1,27 @@
+Copyright (c) 2015 The Go Authors. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+ * Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+ * Neither the name of Google Inc. nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/vendor/github.com/aclements/go-moremath/mathx/beta.go b/vendor/github.com/aclements/go-moremath/mathx/beta.go
new file mode 100644
index 0000000..49f8722
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/beta.go
@@ -0,0 +1,93 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathx
+
+import "math"
+
+func lgamma(x float64) float64 {
+ y, _ := math.Lgamma(x)
+ return y
+}
+
+// Beta returns the value of the complete beta function B(a, b).
+func Beta(a, b float64) float64 {
+ // B(x,y) = Γ(x)Γ(y) / Γ(x+y)
+ return math.Exp(lgamma(a) + lgamma(b) - lgamma(a+b))
+}
+
+// BetaInc returns the value of the regularized incomplete beta
+// function Iₓ(a, b).
+//
+// This is not to be confused with the "incomplete beta function",
+// which can be computed as BetaInc(x, a, b)*Beta(a, b).
+//
+// If x < 0 or x > 1, returns NaN.
+func BetaInc(x, a, b float64) float64 {
+ // Based on Numerical Recipes in C, section 6.4. This uses the
+ // continued fraction definition of I:
+ //
+ // (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...))))))
+ //
+ // where B(a,b) is the beta function and
+ //
+ // d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1))
+ // d_{2m} = m(b-m)x/((a+2m-1)(a+2m))
+ if x < 0 || x > 1 {
+ return math.NaN()
+ }
+ bt := 0.0
+ if 0 < x && x < 1 {
+ // Compute the coefficient before the continued
+ // fraction.
+ bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) +
+ a*math.Log(x) + b*math.Log(1-x))
+ }
+ if x < (a+1)/(a+b+2) {
+ // Compute continued fraction directly.
+ return bt * betacf(x, a, b) / a
+ } else {
+ // Compute continued fraction after symmetry transform.
+ return 1 - bt*betacf(1-x, b, a)/b
+ }
+}
+
+// betacf is the continued fraction component of the regularized
+// incomplete beta function Iₓ(a, b).
+func betacf(x, a, b float64) float64 {
+ const maxIterations = 200
+ const epsilon = 3e-14
+
+ raiseZero := func(z float64) float64 {
+ if math.Abs(z) < math.SmallestNonzeroFloat64 {
+ return math.SmallestNonzeroFloat64
+ }
+ return z
+ }
+
+ c := 1.0
+ d := 1 / raiseZero(1-(a+b)*x/(a+1))
+ h := d
+ for m := 1; m <= maxIterations; m++ {
+ mf := float64(m)
+
+ // Even step of the recurrence.
+ numer := mf * (b - mf) * x / ((a + 2*mf - 1) * (a + 2*mf))
+ d = 1 / raiseZero(1+numer*d)
+ c = raiseZero(1 + numer/c)
+ h *= d * c
+
+ // Odd step of the recurrence.
+ numer = -(a + mf) * (a + b + mf) * x / ((a + 2*mf) * (a + 2*mf + 1))
+ d = 1 / raiseZero(1+numer*d)
+ c = raiseZero(1 + numer/c)
+ hfac := d * c
+ h *= hfac
+
+ if math.Abs(hfac-1) < epsilon {
+ return h
+ }
+ }
+ panic("betainc: a or b too big; failed to converge")
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/choose.go b/vendor/github.com/aclements/go-moremath/mathx/choose.go
new file mode 100644
index 0000000..54dc27c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/choose.go
@@ -0,0 +1,62 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathx
+
+import "math"
+
+const smallFactLimit = 20 // 20! => 62 bits
+var smallFact [smallFactLimit + 1]int64
+
+func init() {
+ smallFact[0] = 1
+ fact := int64(1)
+ for n := int64(1); n <= smallFactLimit; n++ {
+ fact *= n
+ smallFact[n] = fact
+ }
+}
+
+// Choose returns the binomial coefficient of n and k.
+func Choose(n, k int) float64 {
+ if k == 0 || k == n {
+ return 1
+ }
+ if k < 0 || n < k {
+ return 0
+ }
+ if n <= smallFactLimit { // Implies k <= smallFactLimit
+ // It's faster to do several integer multiplications
+ // than it is to do an extra integer division.
+ // Remarkably, this is also faster than pre-computing
+ // Pascal's triangle (presumably because this is very
+ // cache efficient).
+ numer := int64(1)
+ for n1 := int64(n - (k - 1)); n1 <= int64(n); n1++ {
+ numer *= n1
+ }
+ denom := smallFact[k]
+ return float64(numer / denom)
+ }
+
+ return math.Exp(lchoose(n, k))
+}
+
+// Lchoose returns math.Log(Choose(n, k)).
+func Lchoose(n, k int) float64 {
+ if k == 0 || k == n {
+ return 0
+ }
+ if k < 0 || n < k {
+ return math.NaN()
+ }
+ return lchoose(n, k)
+}
+
+func lchoose(n, k int) float64 {
+ a, _ := math.Lgamma(float64(n + 1))
+ b, _ := math.Lgamma(float64(k + 1))
+ c, _ := math.Lgamma(float64(n - k + 1))
+ return a - b - c
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/gamma.go b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
new file mode 100644
index 0000000..d11096e
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
@@ -0,0 +1,96 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathx
+
+import "math"
+
+// GammaInc returns the value of the incomplete gamma function (also
+// known as the regularized gamma function):
+//
+// P(a, x) = 1 / Γ(a) * ∫₀ˣ exp(-t) t**(a-1) dt
+func GammaInc(a, x float64) float64 {
+ // Based on Numerical Recipes in C, section 6.2.
+
+ if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) {
+ return math.NaN()
+ }
+
+ if x < a+1 {
+ // Use the series representation, which converges more
+ // rapidly in this range.
+ return gammaIncSeries(a, x)
+ } else {
+ // Use the continued fraction representation.
+ return 1 - gammaIncCF(a, x)
+ }
+}
+
+// GammaIncComp returns the complement of the incomplete gamma
+// function 1 - GammaInc(a, x). This is more numerically stable for
+// values near 0.
+func GammaIncComp(a, x float64) float64 {
+ if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) {
+ return math.NaN()
+ }
+
+ if x < a+1 {
+ return 1 - gammaIncSeries(a, x)
+ } else {
+ return gammaIncCF(a, x)
+ }
+}
+
+func gammaIncSeries(a, x float64) float64 {
+ const maxIterations = 200
+ const epsilon = 3e-14
+
+ if x == 0 {
+ return 0
+ }
+
+ ap := a
+ del := 1 / a
+ sum := del
+ for n := 0; n < maxIterations; n++ {
+ ap++
+ del *= x / ap
+ sum += del
+ if math.Abs(del) < math.Abs(sum)*epsilon {
+ return sum * math.Exp(-x+a*math.Log(x)-lgamma(a))
+ }
+ }
+ panic("a too large; failed to converge")
+}
+
+func gammaIncCF(a, x float64) float64 {
+ const maxIterations = 200
+ const epsilon = 3e-14
+
+ raiseZero := func(z float64) float64 {
+ if math.Abs(z) < math.SmallestNonzeroFloat64 {
+ return math.SmallestNonzeroFloat64
+ }
+ return z
+ }
+
+ b := x + 1 - a
+ c := math.MaxFloat64
+ d := 1 / b
+ h := d
+
+ for i := 1; i <= maxIterations; i++ {
+ an := -float64(i) * (float64(i) - a)
+ b += 2
+ d = raiseZero(an*d + b)
+ c = raiseZero(b + an/c)
+ d = 1 / d
+ del := d * c
+ h *= del
+ if math.Abs(del-1) < epsilon {
+ return math.Exp(-x+a*math.Log(x)-lgamma(a)) * h
+ }
+ }
+ panic("a too large; failed to converge")
+}
diff --git a/vendor/github.com/aclements/go-moremath/mathx/package.go b/vendor/github.com/aclements/go-moremath/mathx/package.go
new file mode 100644
index 0000000..9d5de0d
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/package.go
@@ -0,0 +1,11 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package mathx implements special functions not provided by the
+// standard math package.
+package mathx // import "github.com/aclements/go-moremath/mathx"
+
+import "math"
+
+var nan = math.NaN()
diff --git a/vendor/github.com/aclements/go-moremath/mathx/sign.go b/vendor/github.com/aclements/go-moremath/mathx/sign.go
new file mode 100644
index 0000000..372e92f
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/sign.go
@@ -0,0 +1,18 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathx
+
+// Sign returns the sign of x: -1 if x < 0, 0 if x == 0, 1 if x > 0.
+// If x is NaN, it returns NaN.
+func Sign(x float64) float64 {
+ if x == 0 {
+ return 0
+ } else if x < 0 {
+ return -1
+ } else if x > 0 {
+ return 1
+ }
+ return nan
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/alg.go b/vendor/github.com/aclements/go-moremath/stats/alg.go
new file mode 100644
index 0000000..f704cef
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/alg.go
@@ -0,0 +1,114 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+// Miscellaneous helper algorithms
+
+import (
+ "fmt"
+
+ "github.com/aclements/go-moremath/mathx"
+)
+
+func maxint(a, b int) int {
+ if a > b {
+ return a
+ }
+ return b
+}
+
+func minint(a, b int) int {
+ if a < b {
+ return a
+ }
+ return b
+}
+
+func sumint(xs []int) int {
+ sum := 0
+ for _, x := range xs {
+ sum += x
+ }
+ return sum
+}
+
+// bisect returns an x in [low, high] such that |f(x)| <= tolerance
+// using the bisection method.
+//
+// f(low) and f(high) must have opposite signs.
+//
+// If f does not have a root in this interval (e.g., it is
+// discontiguous), this returns the X of the apparent discontinuity
+// and false.
+func bisect(f func(float64) float64, low, high, tolerance float64) (float64, bool) {
+ flow, fhigh := f(low), f(high)
+ if -tolerance <= flow && flow <= tolerance {
+ return low, true
+ }
+ if -tolerance <= fhigh && fhigh <= tolerance {
+ return high, true
+ }
+ if mathx.Sign(flow) == mathx.Sign(fhigh) {
+ panic(fmt.Sprintf("root of f is not bracketed by [low, high]; f(%g)=%g f(%g)=%g", low, flow, high, fhigh))
+ }
+ for {
+ mid := (high + low) / 2
+ fmid := f(mid)
+ if -tolerance <= fmid && fmid <= tolerance {
+ return mid, true
+ }
+ if mid == high || mid == low {
+ return mid, false
+ }
+ if mathx.Sign(fmid) == mathx.Sign(flow) {
+ low = mid
+ flow = fmid
+ } else {
+ high = mid
+ fhigh = fmid
+ }
+ }
+}
+
+// bisectBool implements the bisection method on a boolean function.
+// It returns x1, x2 ∈ [low, high], x1 < x2 such that f(x1) != f(x2)
+// and x2 - x1 <= xtol.
+//
+// If f(low) == f(high), it panics.
+func bisectBool(f func(float64) bool, low, high, xtol float64) (x1, x2 float64) {
+ flow, fhigh := f(low), f(high)
+ if flow == fhigh {
+ panic(fmt.Sprintf("root of f is not bracketed by [low, high]; f(%g)=%v f(%g)=%v", low, flow, high, fhigh))
+ }
+ for {
+ if high-low <= xtol {
+ return low, high
+ }
+ mid := (high + low) / 2
+ if mid == high || mid == low {
+ return low, high
+ }
+ fmid := f(mid)
+ if fmid == flow {
+ low = mid
+ flow = fmid
+ } else {
+ high = mid
+ fhigh = fmid
+ }
+ }
+}
+
+// series returns the sum of the series f(0), f(1), ...
+//
+// This implementation is fast, but subject to round-off error.
+func series(f func(float64) float64) float64 {
+ y, yp := 0.0, 1.0
+ for n := 0.0; y != yp; n++ {
+ yp = y
+ y += f(n)
+ }
+ return y
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/deltadist.go b/vendor/github.com/aclements/go-moremath/stats/deltadist.go
new file mode 100644
index 0000000..bb3ba3f
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/deltadist.go
@@ -0,0 +1,57 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+// DeltaDist is the Dirac delta function, centered at T, with total
+// area 1.
+//
+// The CDF of the Dirac delta function is the Heaviside step function,
+// centered at T. Specifically, f(T) == 1.
+type DeltaDist struct {
+ T float64
+}
+
+func (d DeltaDist) PDF(x float64) float64 {
+ if x == d.T {
+ return inf
+ }
+ return 0
+}
+
+func (d DeltaDist) pdfEach(xs []float64) []float64 {
+ res := make([]float64, len(xs))
+ for i, x := range xs {
+ if x == d.T {
+ res[i] = inf
+ }
+ }
+ return res
+}
+
+func (d DeltaDist) CDF(x float64) float64 {
+ if x >= d.T {
+ return 1
+ }
+ return 0
+}
+
+func (d DeltaDist) cdfEach(xs []float64) []float64 {
+ res := make([]float64, len(xs))
+ for i, x := range xs {
+ res[i] = d.CDF(x)
+ }
+ return res
+}
+
+func (d DeltaDist) InvCDF(y float64) float64 {
+ if y < 0 || y > 1 {
+ return nan
+ }
+ return d.T
+}
+
+func (d DeltaDist) Bounds() (float64, float64) {
+ return d.T - 1, d.T + 1
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/dist.go b/vendor/github.com/aclements/go-moremath/stats/dist.go
new file mode 100644
index 0000000..048477d
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/dist.go
@@ -0,0 +1,210 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import "math/rand"
+
+// A DistCommon is a statistical distribution. DistCommon is a base
+// interface provided by both continuous and discrete distributions.
+type DistCommon interface {
+ // CDF returns the cumulative probability Pr[X <= x].
+ //
+ // For continuous distributions, the CDF is the integral of
+ // the PDF from -inf to x.
+ //
+ // For discrete distributions, the CDF is the sum of the PMF
+ // at all defined points from -inf to x, inclusive. Note that
+ // the CDF of a discrete distribution is defined for the whole
+ // real line (unlike the PMF) but has discontinuities where
+ // the PMF is non-zero.
+ //
+ // The CDF is a monotonically increasing function and has a
+ // domain of all real numbers. If the distribution has bounded
+ // support, it has a range of [0, 1]; otherwise it has a range
+ // of (0, 1). Finally, CDF(-inf)==0 and CDF(inf)==1.
+ CDF(x float64) float64
+
+ // Bounds returns reasonable bounds for this distribution's
+ // PDF/PMF and CDF. The total weight outside of these bounds
+ // should be approximately 0.
+ //
+ // For a discrete distribution, both bounds are integer
+ // multiples of Step().
+ //
+ // If this distribution has finite support, it returns exact
+ // bounds l, h such that CDF(l')=0 for all l' < l and
+ // CDF(h')=1 for all h' >= h.
+ Bounds() (float64, float64)
+}
+
+// A Dist is a continuous statistical distribution.
+type Dist interface {
+ DistCommon
+
+ // PDF returns the value of the probability density function
+ // of this distribution at x.
+ PDF(x float64) float64
+}
+
+// A DiscreteDist is a discrete statistical distribution.
+//
+// Most discrete distributions are defined only at integral values of
+// the random variable. However, some are defined at other intervals,
+// so this interface takes a float64 value for the random variable.
+// The probability mass function rounds down to the nearest defined
+// point. Note that float64 values can exactly represent integer
+// values between ±2**53, so this generally shouldn't be an issue for
+// integer-valued distributions (likewise, for half-integer-valued
+// distributions, float64 can exactly represent all values between
+// ±2**52).
+type DiscreteDist interface {
+ DistCommon
+
+ // PMF returns the value of the probability mass function
+ // Pr[X = x'], where x' is x rounded down to the nearest
+ // defined point on the distribution.
+ //
+ // Note for implementers: for integer-valued distributions,
+ // round x using int(math.Floor(x)). Do not use int(x), since
+ // that truncates toward zero (unless all x <= 0 are handled
+ // the same).
+ PMF(x float64) float64
+
+ // Step returns s, where the distribution is defined for sℕ.
+ Step() float64
+}
+
+// TODO: Add a Support method for finite support distributions? Or
+// maybe just another return value from Bounds indicating that the
+// bounds are exact?
+
+// TODO: Plot method to return a pre-configured Plot object with
+// reasonable bounds and an integral function? Have to distinguish
+// PDF/CDF/InvCDF. Three methods? Argument?
+//
+// Doesn't have to be a method of Dist. Could be just a function that
+// takes a Dist and uses Bounds.
+
+// InvCDF returns the inverse CDF function of the given distribution
+// (also known as the quantile function or the percent point
+// function). This is a function f such that f(dist.CDF(x)) == x. If
+// dist.CDF is only weakly monotonic (that it, there are intervals
+// over which it is constant) and y > 0, f returns the smallest x that
+// satisfies this condition. In general, the inverse CDF is not
+// well-defined for y==0, but for convenience if y==0, f returns the
+// largest x that satisfies this condition. For distributions with
+// infinite support both the largest and smallest x are -Inf; however,
+// for distributions with finite support, this is the lower bound of
+// the support.
+//
+// If y < 0 or y > 1, f returns NaN.
+//
+// If dist implements InvCDF(float64) float64, this returns that
+// method. Otherwise, it returns a function that uses a generic
+// numerical method to construct the inverse CDF at y by finding x
+// such that dist.CDF(x) == y. This may have poor precision around
+// points of discontinuity, including f(0) and f(1).
+func InvCDF(dist DistCommon) func(y float64) (x float64) {
+ type invCDF interface {
+ InvCDF(float64) float64
+ }
+ if dist, ok := dist.(invCDF); ok {
+ return dist.InvCDF
+ }
+
+ // Otherwise, use a numerical algorithm.
+ //
+ // TODO: For discrete distributions, use the step size to
+ // inform this computation.
+ return func(y float64) (x float64) {
+ const almostInf = 1e100
+ const xtol = 1e-16
+
+ if y < 0 || y > 1 {
+ return nan
+ } else if y == 0 {
+ l, _ := dist.Bounds()
+ if dist.CDF(l) == 0 {
+ // Finite support
+ return l
+ } else {
+ // Infinite support
+ return -inf
+ }
+ } else if y == 1 {
+ _, h := dist.Bounds()
+ if dist.CDF(h) == 1 {
+ // Finite support
+ return h
+ } else {
+ // Infinite support
+ return inf
+ }
+ }
+
+ // Find loX, hiX for which cdf(loX) < y <= cdf(hiX).
+ var loX, loY, hiX, hiY float64
+ x1, y1 := 0.0, dist.CDF(0)
+ xdelta := 1.0
+ if y1 < y {
+ hiX, hiY = x1, y1
+ for hiY < y && hiX != inf {
+ loX, loY, hiX = hiX, hiY, hiX+xdelta
+ hiY = dist.CDF(hiX)
+ xdelta *= 2
+ }
+ } else {
+ loX, loY = x1, y1
+ for y <= loY && loX != -inf {
+ hiX, hiY, loX = loX, loY, loX-xdelta
+ loY = dist.CDF(loX)
+ xdelta *= 2
+ }
+ }
+ if loX == -inf {
+ return loX
+ } else if hiX == inf {
+ return hiX
+ }
+
+ // Use bisection on the interval to find the smallest
+ // x at which cdf(x) <= y.
+ _, x = bisectBool(func(x float64) bool {
+ return dist.CDF(x) < y
+ }, loX, hiX, xtol)
+ return
+ }
+}
+
+// Rand returns a random number generator that draws from the given
+// distribution. The returned generator takes an optional source of
+// randomness; if this is nil, it uses the default global source.
+//
+// If dist implements Rand(*rand.Rand) float64, Rand returns that
+// method. Otherwise, it returns a generic generator based on dist's
+// inverse CDF (which may in turn use an efficient implementation or a
+// generic numerical implementation; see InvCDF).
+func Rand(dist DistCommon) func(*rand.Rand) float64 {
+ type distRand interface {
+ Rand(*rand.Rand) float64
+ }
+ if dist, ok := dist.(distRand); ok {
+ return dist.Rand
+ }
+
+ // Otherwise, use a generic algorithm.
+ inv := InvCDF(dist)
+ return func(r *rand.Rand) float64 {
+ var y float64
+ for y == 0 {
+ if r == nil {
+ y = rand.Float64()
+ } else {
+ y = r.Float64()
+ }
+ }
+ return inv(y)
+ }
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/hist.go b/vendor/github.com/aclements/go-moremath/stats/hist.go
new file mode 100644
index 0000000..8578c07
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/hist.go
@@ -0,0 +1,66 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import "math"
+
+// TODO: Implement histograms on top of scales.
+
+type Histogram interface {
+ // Add adds a sample with value x to histogram h.
+ Add(x float64)
+
+ // Counts returns the number of samples less than the lowest
+ // bin, a slice of the number of samples in each bin,
+ // and the number of samples greater than the highest bin.
+ Counts() (under uint, counts []uint, over uint)
+
+ // BinToValue returns the value that would appear at the given
+ // bin index.
+ //
+ // For integral values of bin, BinToValue returns the lower
+ // bound of bin. That is, a sample value x will be in bin if
+ // bin is integral and
+ //
+ // BinToValue(bin) <= x < BinToValue(bin + 1)
+ //
+ // For non-integral values of bin, BinToValue interpolates
+ // between the lower and upper bounds of math.Floor(bin).
+ //
+ // BinToValue is undefined if bin > 1 + the number of bins.
+ BinToValue(bin float64) float64
+}
+
+// HistogramQuantile returns the x such that n*q samples in hist are
+// <= x, assuming values are distibuted within each bin according to
+// hist's distribution.
+//
+// If the q'th sample falls below the lowest bin or above the highest
+// bin, returns NaN.
+func HistogramQuantile(hist Histogram, q float64) float64 {
+ under, counts, over := hist.Counts()
+ total := under + over
+ for _, count := range counts {
+ total += count
+ }
+
+ goal := uint(float64(total) * q)
+ if goal <= under || goal > total-over {
+ return math.NaN()
+ }
+ for bin, count := range counts {
+ if count > goal {
+ return hist.BinToValue(float64(bin) + float64(goal)/float64(count))
+ }
+ goal -= count
+ }
+ panic("goal count not reached")
+}
+
+// HistogramIQR returns the interquartile range of the samples in
+// hist.
+func HistogramIQR(hist Histogram) float64 {
+ return HistogramQuantile(hist, 0.75) - HistogramQuantile(hist, 0.25)
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/hypergdist.go b/vendor/github.com/aclements/go-moremath/stats/hypergdist.go
new file mode 100644
index 0000000..1ea05e2
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/hypergdist.go
@@ -0,0 +1,101 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+
+ "github.com/aclements/go-moremath/mathx"
+)
+
+// HypergeometicDist is a hypergeometric distribution.
+type HypergeometicDist struct {
+ // N is the size of the population. N >= 0.
+ N int
+
+ // K is the number of successes in the population. 0 <= K <= N.
+ K int
+
+ // Draws is the number of draws from the population. This is
+ // usually written "n", but is called Draws here because of
+ // limitations on Go identifier naming. 0 <= Draws <= N.
+ Draws int
+}
+
+// PMF is the probability of getting exactly int(k) successes in
+// d.Draws draws with replacement from a population of size d.N that
+// contains exactly d.K successes.
+func (d HypergeometicDist) PMF(k float64) float64 {
+ ki := int(math.Floor(k))
+ l, h := d.bounds()
+ if ki < l || ki > h {
+ return 0
+ }
+ return d.pmf(ki)
+}
+
+func (d HypergeometicDist) pmf(k int) float64 {
+ return math.Exp(mathx.Lchoose(d.K, k) + mathx.Lchoose(d.N-d.K, d.Draws-k) - mathx.Lchoose(d.N, d.Draws))
+}
+
+// CDF is the probability of getting int(k) or fewer successes in
+// d.Draws draws with replacement from a population of size d.N that
+// contains exactly d.K successes.
+func (d HypergeometicDist) CDF(k float64) float64 {
+ // Based on Klotz, A Computational Approach to Statistics.
+ ki := int(math.Floor(k))
+ l, h := d.bounds()
+ if ki < l {
+ return 0
+ } else if ki >= h {
+ return 1
+ }
+ // Use symmetry to compute the smaller sum.
+ flip := false
+ if ki > (d.Draws+1)/(d.N+1)*(d.K+1) {
+ flip = true
+ ki = d.K - ki - 1
+ d.Draws = d.N - d.Draws
+ }
+ p := d.pmf(ki) * d.sum(ki)
+ if flip {
+ p = 1 - p
+ }
+ return p
+}
+
+func (d HypergeometicDist) sum(k int) float64 {
+ const epsilon = 1e-14
+ sum, ak := 1.0, 1.0
+ L := maxint(0, d.Draws+d.K-d.N)
+ for dk := 1; dk <= k-L && ak/sum > epsilon; dk++ {
+ ak *= float64(1+k-dk) / float64(d.Draws-k+dk)
+ ak *= float64(d.N-d.K-d.Draws+k+1-dk) / float64(d.K-k+dk)
+ sum += ak
+ }
+ return sum
+}
+
+func (d HypergeometicDist) bounds() (int, int) {
+ return maxint(0, d.Draws+d.K-d.N), minint(d.Draws, d.K)
+}
+
+func (d HypergeometicDist) Bounds() (float64, float64) {
+ l, h := d.bounds()
+ return float64(l), float64(h)
+}
+
+func (d HypergeometicDist) Step() float64 {
+ return 1
+}
+
+func (d HypergeometicDist) Mean() float64 {
+ return float64(d.Draws*d.K) / float64(d.N)
+}
+
+func (d HypergeometicDist) Variance() float64 {
+ return float64(d.Draws*d.K*(d.N-d.K)*(d.N-d.Draws)) /
+ float64(d.N*d.N*(d.N-1))
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kde.go b/vendor/github.com/aclements/go-moremath/stats/kde.go
new file mode 100644
index 0000000..6bcbd49
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kde.go
@@ -0,0 +1,350 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "fmt"
+ "math"
+)
+
+// A KDE is a distribution that estimates the underlying distribution
+// of a Sample using kernel density estimation.
+//
+// Kernel density estimation is a method for constructing an estimate
+// ƒ̂(x) of a unknown distribution ƒ(x) given a sample from that
+// distribution. Unlike many techniques, kernel density estimation is
+// non-parametric: in general, it doesn't assume any particular true
+// distribution (note, however, that the resulting distribution
+// depends deeply on the selected bandwidth, and many bandwidth
+// estimation techniques assume normal reference rules).
+//
+// A kernel density estimate is similar to a histogram, except that it
+// is a smooth probability estimate and does not require choosing a
+// bin size and discretizing the data.
+//
+// Sample is the only required field. All others have reasonable
+// defaults.
+type KDE struct {
+ // Sample is the data sample underlying this KDE.
+ Sample Sample
+
+ // Kernel is the kernel to use for the KDE.
+ Kernel KDEKernel
+
+ // Bandwidth is the bandwidth to use for the KDE.
+ //
+ // If this is zero, the bandwidth is computed from the
+ // provided data using a default bandwidth estimator
+ // (currently BandwidthScott).
+ Bandwidth float64
+
+ // BoundaryMethod is the boundary correction method to use for
+ // the KDE. The default value is BoundaryReflect; however, the
+ // default bounds are effectively +/-inf, which is equivalent
+ // to performing no boundary correction.
+ BoundaryMethod KDEBoundaryMethod
+
+ // [BoundaryMin, BoundaryMax) specify a bounded support for
+ // the KDE. If both are 0 (their default values), they are
+ // treated as +/-inf.
+ //
+ // To specify a half-bounded support, set Min to math.Inf(-1)
+ // or Max to math.Inf(1).
+ BoundaryMin float64
+ BoundaryMax float64
+}
+
+// BandwidthSilverman is a bandwidth estimator implementing
+// Silverman's Rule of Thumb. It's fast, but not very robust to
+// outliers as it assumes data is approximately normal.
+//
+// Silverman, B. W. (1986) Density Estimation.
+func BandwidthSilverman(data interface {
+ StdDev() float64
+ Weight() float64
+}) float64 {
+ return 1.06 * data.StdDev() * math.Pow(data.Weight(), -1.0/5)
+}
+
+// BandwidthScott is a bandwidth estimator implementing Scott's Rule.
+// This is generally robust to outliers: it chooses the minimum
+// between the sample's standard deviation and an robust estimator of
+// a Gaussian distribution's standard deviation.
+//
+// Scott, D. W. (1992) Multivariate Density Estimation: Theory,
+// Practice, and Visualization.
+func BandwidthScott(data interface {
+ StdDev() float64
+ Weight() float64
+ Quantile(float64) float64
+}) float64 {
+ iqr := data.Quantile(0.75) - data.Quantile(0.25)
+ hScale := 1.06 * math.Pow(data.Weight(), -1.0/5)
+ stdDev := data.StdDev()
+ if stdDev < iqr/1.349 {
+ // Use Silverman's Rule of Thumb
+ return hScale * stdDev
+ } else {
+ // Use IQR/1.349 as a robust estimator of the standard
+ // deviation of a Gaussian distribution.
+ return hScale * (iqr / 1.349)
+ }
+}
+
+// TODO(austin) Implement bandwidth estimator from Botev, Grotowski,
+// Kroese. (2010) Kernel Density Estimation via Diffusion.
+
+// KDEKernel represents a kernel to use for a KDE.
+type KDEKernel int
+
+//go:generate stringer -type=KDEKernel
+
+const (
+ // An EpanechnikovKernel is a smooth kernel with bounded
+ // support. As a result, the KDE will also have bounded
+ // support. It is "optimal" in the sense that it minimizes the
+ // asymptotic mean integrated squared error (AMISE).
+ EpanechnikovKernel KDEKernel = iota
+
+ // A GaussianKernel is a Gaussian (normal) kernel.
+ GaussianKernel
+
+ // A DeltaKernel is a Dirac delta function. The PDF of such a
+ // KDE is not well-defined, but the CDF will represent each
+ // sample as an instantaneous increase. This kernel ignores
+ // bandwidth and never requires boundary correction.
+ DeltaKernel
+)
+
+// KDEBoundaryMethod represents a boundary correction method for
+// constructing a KDE with bounded support.
+type KDEBoundaryMethod int
+
+//go:generate stringer -type=KDEBoundaryMethod
+
+const (
+ // BoundaryReflect reflects the density estimate at the
+ // boundaries. For example, for a KDE with support [0, inf),
+ // this is equivalent to ƒ̂ᵣ(x)=ƒ̂(x)+ƒ̂(-x) for x>=0. This is a
+ // simple and fast technique, but enforces that ƒ̂ᵣ'(0)=0, so
+ // it may not be applicable to all distributions.
+ BoundaryReflect KDEBoundaryMethod = iota
+)
+
+type kdeKernel interface {
+ pdfEach(xs []float64) []float64
+ cdfEach(xs []float64) []float64
+}
+
+func (k *KDE) prepare() (kdeKernel, bool) {
+ // Compute bandwidth.
+ if k.Bandwidth == 0 {
+ k.Bandwidth = BandwidthScott(k.Sample)
+ }
+
+ // Construct kernel.
+ kernel := kdeKernel(nil)
+ switch k.Kernel {
+ default:
+ panic(fmt.Sprint("unknown kernel", k))
+ case EpanechnikovKernel:
+ kernel = epanechnikovKernel{k.Bandwidth}
+ case GaussianKernel:
+ kernel = NormalDist{0, k.Bandwidth}
+ case DeltaKernel:
+ kernel = DeltaDist{0}
+ }
+
+ // Use boundary correction?
+ bc := k.BoundaryMin != 0 || k.BoundaryMax != 0
+
+ return kernel, bc
+}
+
+// TODO: For KDEs of histograms, make histograms able to create a
+// weighted Sample and simply require the caller to provide a
+// good bandwidth from a StreamStats.
+
+// normalizedXs returns x - kde.Sample.Xs. Evaluating kernels shifted
+// by kde.Sample.Xs all at x is equivalent to evaluating one unshifted
+// kernel at x - kde.Sample.Xs.
+func (kde *KDE) normalizedXs(x float64) []float64 {
+ txs := make([]float64, len(kde.Sample.Xs))
+ for i, xi := range kde.Sample.Xs {
+ txs[i] = x - xi
+ }
+ return txs
+}
+
+func (kde *KDE) PDF(x float64) float64 {
+ kernel, bc := kde.prepare()
+
+ // Apply boundary
+ if bc && (x < kde.BoundaryMin || x >= kde.BoundaryMax) {
+ return 0
+ }
+
+ y := func(x float64) float64 {
+ // Shift kernel to each of kde.xs and evaluate at x
+ ys := kernel.pdfEach(kde.normalizedXs(x))
+
+ // Kernel samples are weighted according to the weights of xs
+ wys := Sample{Xs: ys, Weights: kde.Sample.Weights}
+
+ return wys.Sum() / wys.Weight()
+ }
+ if !bc {
+ return y(x)
+ }
+ switch kde.BoundaryMethod {
+ default:
+ panic("unknown boundary correction method")
+ case BoundaryReflect:
+ if math.IsInf(kde.BoundaryMax, 1) {
+ return y(x) + y(2*kde.BoundaryMin-x)
+ } else if math.IsInf(kde.BoundaryMin, -1) {
+ return y(x) + y(2*kde.BoundaryMax-x)
+ } else {
+ d := 2 * (kde.BoundaryMax - kde.BoundaryMin)
+ w := 2 * (x - kde.BoundaryMin)
+ return series(func(n float64) float64 {
+ // Points >= x
+ return y(x+n*d) + y(x+n*d-w)
+ }) + series(func(n float64) float64 {
+ // Points < x
+ return y(x-(n+1)*d+w) + y(x-(n+1)*d)
+ })
+ }
+ }
+}
+
+func (kde *KDE) CDF(x float64) float64 {
+ kernel, bc := kde.prepare()
+
+ // Apply boundary
+ if bc {
+ if x < kde.BoundaryMin {
+ return 0
+ } else if x >= kde.BoundaryMax {
+ return 1
+ }
+ }
+
+ y := func(x float64) float64 {
+ // Shift kernel integral to each of cdf.xs and evaluate at x
+ ys := kernel.cdfEach(kde.normalizedXs(x))
+
+ // Kernel samples are weighted according to the weights of xs
+ wys := Sample{Xs: ys, Weights: kde.Sample.Weights}
+
+ return wys.Sum() / wys.Weight()
+ }
+ if !bc {
+ return y(x)
+ }
+ switch kde.BoundaryMethod {
+ default:
+ panic("unknown boundary correction method")
+ case BoundaryReflect:
+ if math.IsInf(kde.BoundaryMax, 1) {
+ return y(x) - y(2*kde.BoundaryMin-x)
+ } else if math.IsInf(kde.BoundaryMin, -1) {
+ return y(x) + (1 - y(2*kde.BoundaryMax-x))
+ } else {
+ d := 2 * (kde.BoundaryMax - kde.BoundaryMin)
+ w := 2 * (x - kde.BoundaryMin)
+ return series(func(n float64) float64 {
+ // Windows >= x-w
+ return y(x+n*d) - y(x+n*d-w)
+ }) + series(func(n float64) float64 {
+ // Windows < x-w
+ return y(x-(n+1)*d) - y(x-(n+1)*d-w)
+ })
+ }
+ }
+}
+
+func (kde *KDE) Bounds() (low float64, high float64) {
+ _, bc := kde.prepare()
+
+ // TODO(austin) If this KDE came from a histogram, we'd better
+ // not sample at a significantly higher rate than the
+ // histogram. Maybe we want to just return the bounds of the
+ // histogram?
+
+ // TODO(austin) It would be nice if this could be instructed
+ // to include all original data points, even if they are in
+ // the tail. Probably that should just be up to the caller to
+ // pass an axis derived from the bounds of the original data.
+
+ // Use the lowest and highest samples as starting points
+ lowX, highX := kde.Sample.Bounds()
+ if lowX == highX {
+ lowX -= 1
+ highX += 1
+ }
+
+ // Find the end points that contain 99% of the CDF's weight.
+ // Since bisect requires that the root be bracketed, start by
+ // expanding our range if necessary. TODO(austin) This can
+ // definitely be done faster.
+ const (
+ lowY = 0.005
+ highY = 0.995
+ tolerance = 0.001
+ )
+ for kde.CDF(lowX) > lowY {
+ lowX -= highX - lowX
+ }
+ for kde.CDF(highX) < highY {
+ highX += highX - lowX
+ }
+ // Explicitly accept discontinuities, since we may be using a
+ // discontiguous kernel.
+ low, _ = bisect(func(x float64) float64 { return kde.CDF(x) - lowY }, lowX, highX, tolerance)
+ high, _ = bisect(func(x float64) float64 { return kde.CDF(x) - highY }, lowX, highX, tolerance)
+
+ // Expand width by 20% to give some margins
+ width := high - low
+ low, high = low-0.1*width, high+0.1*width
+
+ // Limit to bounds
+ if bc {
+ low = math.Max(low, kde.BoundaryMin)
+ high = math.Min(high, kde.BoundaryMax)
+ }
+
+ return
+}
+
+type epanechnikovKernel struct {
+ h float64
+}
+
+func (d epanechnikovKernel) pdfEach(xs []float64) []float64 {
+ ys := make([]float64, len(xs))
+ a := 0.75 / d.h
+ invhh := 1 / (d.h * d.h)
+ for i, x := range xs {
+ if -d.h < x && x < d.h {
+ ys[i] = a * (1 - x*x*invhh)
+ }
+ }
+ return ys
+}
+
+func (d epanechnikovKernel) cdfEach(xs []float64) []float64 {
+ ys := make([]float64, len(xs))
+ invh := 1 / d.h
+ for i, x := range xs {
+ if x > d.h {
+ ys[i] = 1
+ } else if x > -d.h {
+ u := x * invh
+ ys[i] = 0.25 * (2 + 3*u - u*u*u)
+ }
+ }
+ return ys
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go b/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go
new file mode 100644
index 0000000..e01d0e7
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kdeboundarymethod_string.go
@@ -0,0 +1,16 @@
+// generated by stringer -type=KDEBoundaryMethod; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _KDEBoundaryMethod_name = "BoundaryReflect"
+
+var _KDEBoundaryMethod_index = [...]uint8{0, 15}
+
+func (i KDEBoundaryMethod) String() string {
+ if i < 0 || i+1 >= KDEBoundaryMethod(len(_KDEBoundaryMethod_index)) {
+ return fmt.Sprintf("KDEBoundaryMethod(%d)", i)
+ }
+ return _KDEBoundaryMethod_name[_KDEBoundaryMethod_index[i]:_KDEBoundaryMethod_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go b/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go
new file mode 100644
index 0000000..20a23fd
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/kdekernel_string.go
@@ -0,0 +1,16 @@
+// generated by stringer -type=KDEKernel; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _KDEKernel_name = "GaussianKernelDeltaKernel"
+
+var _KDEKernel_index = [...]uint8{0, 14, 25}
+
+func (i KDEKernel) String() string {
+ if i < 0 || i+1 >= KDEKernel(len(_KDEKernel_index)) {
+ return fmt.Sprintf("KDEKernel(%d)", i)
+ }
+ return _KDEKernel_name[_KDEKernel_index[i]:_KDEKernel_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/linearhist.go b/vendor/github.com/aclements/go-moremath/stats/linearhist.go
new file mode 100644
index 0000000..c36335a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/linearhist.go
@@ -0,0 +1,42 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+// LinearHist is a Histogram with uniformly-sized bins.
+type LinearHist struct {
+ min, max, delta float64
+ low, high uint
+ bins []uint
+}
+
+// NewLinearHist returns an empty histogram with nbins uniformly-sized
+// bins spanning [min, max].
+func NewLinearHist(min, max float64, nbins int) *LinearHist {
+ delta := float64(nbins) / (max - min)
+ return &LinearHist{min, max, delta, 0, 0, make([]uint, nbins)}
+}
+
+func (h *LinearHist) bin(x float64) int {
+ return int(h.delta * (x - h.min))
+}
+
+func (h *LinearHist) Add(x float64) {
+ bin := h.bin(x)
+ if bin < 0 {
+ h.low++
+ } else if bin >= len(h.bins) {
+ h.high++
+ } else {
+ h.bins[bin]++
+ }
+}
+
+func (h *LinearHist) Counts() (uint, []uint, uint) {
+ return h.low, h.bins, h.high
+}
+
+func (h *LinearHist) BinToValue(bin float64) float64 {
+ return h.min + bin*h.delta
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go b/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go
new file mode 100644
index 0000000..ab0f26c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/locationhypothesis_string.go
@@ -0,0 +1,17 @@
+// generated by stringer -type LocationHypothesis; DO NOT EDIT
+
+package stats
+
+import "fmt"
+
+const _LocationHypothesis_name = "LocationLessLocationDiffersLocationGreater"
+
+var _LocationHypothesis_index = [...]uint8{0, 12, 27, 42}
+
+func (i LocationHypothesis) String() string {
+ i -= -1
+ if i < 0 || i+1 >= LocationHypothesis(len(_LocationHypothesis_index)) {
+ return fmt.Sprintf("LocationHypothesis(%d)", i+-1)
+ }
+ return _LocationHypothesis_name[_LocationHypothesis_index[i]:_LocationHypothesis_index[i+1]]
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/loghist.go b/vendor/github.com/aclements/go-moremath/stats/loghist.go
new file mode 100644
index 0000000..937f62a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/loghist.go
@@ -0,0 +1,83 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import "math"
+
+// LogHist is a Histogram with logarithmically-spaced bins.
+type LogHist struct {
+ b int
+ m float64
+ mOverLogb float64
+ low, high uint
+ bins []uint
+}
+
+// NewLogHist returns an empty logarithmic histogram with bins for
+// integral values of m * log_b(x) up to x = max.
+func NewLogHist(b int, m float64, max float64) *LogHist {
+ // TODO(austin) Minimum value as well? If the samples are
+ // actually integral, having fractional bin boundaries can
+ // mess up smoothing.
+ mOverLogb := m / math.Log(float64(b))
+ nbins := int(math.Ceil(mOverLogb * math.Log(max)))
+ return &LogHist{b: b, m: m, mOverLogb: mOverLogb, low: 0, high: 0, bins: make([]uint, nbins)}
+}
+
+func (h *LogHist) bin(x float64) int {
+ return int(h.mOverLogb * math.Log(x))
+}
+
+func (h *LogHist) Add(x float64) {
+ bin := h.bin(x)
+ if bin < 0 {
+ h.low++
+ } else if bin >= len(h.bins) {
+ h.high++
+ } else {
+ h.bins[bin]++
+ }
+}
+
+func (h *LogHist) Counts() (uint, []uint, uint) {
+ return h.low, h.bins, h.high
+}
+
+func (h *LogHist) BinToValue(bin float64) float64 {
+ return math.Pow(float64(h.b), bin/h.m)
+}
+
+func (h *LogHist) At(x float64) float64 {
+ bin := h.bin(x)
+ if bin < 0 || bin >= len(h.bins) {
+ return 0
+ }
+ return float64(h.bins[bin])
+}
+
+func (h *LogHist) Bounds() (float64, float64) {
+ // XXX Plot will plot this on a linear axis. Maybe this
+ // should be able to return the natural axis?
+ // Maybe then we could also give it the bins for the tics.
+ lowbin := 0
+ if h.low == 0 {
+ for bin, count := range h.bins {
+ if count > 0 {
+ lowbin = bin
+ break
+ }
+ }
+ }
+ highbin := len(h.bins)
+ if h.high == 0 {
+ for bin := range h.bins {
+ if h.bins[len(h.bins)-bin-1] > 0 {
+ highbin = len(h.bins) - bin
+ break
+ }
+ }
+ }
+ return h.BinToValue(float64(lowbin)), h.BinToValue(float64(highbin))
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/normaldist.go b/vendor/github.com/aclements/go-moremath/stats/normaldist.go
new file mode 100644
index 0000000..d00f96a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/normaldist.go
@@ -0,0 +1,141 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+ "math/rand"
+)
+
+// NormalDist is a normal (Gaussian) distribution with mean Mu and
+// standard deviation Sigma.
+type NormalDist struct {
+ Mu, Sigma float64
+}
+
+// StdNormal is the standard normal distribution (Mu = 0, Sigma = 1)
+var StdNormal = NormalDist{0, 1}
+
+// 1/sqrt(2 * pi)
+const invSqrt2Pi = 0.39894228040143267793994605993438186847585863116493465766592583
+
+func (n NormalDist) PDF(x float64) float64 {
+ z := x - n.Mu
+ return math.Exp(-z*z/(2*n.Sigma*n.Sigma)) * invSqrt2Pi / n.Sigma
+}
+
+func (n NormalDist) pdfEach(xs []float64) []float64 {
+ res := make([]float64, len(xs))
+ if n.Mu == 0 && n.Sigma == 1 {
+ // Standard normal fast path
+ for i, x := range xs {
+ res[i] = math.Exp(-x*x/2) * invSqrt2Pi
+ }
+ } else {
+ a := -1 / (2 * n.Sigma * n.Sigma)
+ b := invSqrt2Pi / n.Sigma
+ for i, x := range xs {
+ z := x - n.Mu
+ res[i] = math.Exp(z*z*a) * b
+ }
+ }
+ return res
+}
+
+func (n NormalDist) CDF(x float64) float64 {
+ return math.Erfc(-(x-n.Mu)/(n.Sigma*math.Sqrt2)) / 2
+}
+
+func (n NormalDist) cdfEach(xs []float64) []float64 {
+ res := make([]float64, len(xs))
+ a := 1 / (n.Sigma * math.Sqrt2)
+ for i, x := range xs {
+ res[i] = math.Erfc(-(x-n.Mu)*a) / 2
+ }
+ return res
+}
+
+func (n NormalDist) InvCDF(p float64) (x float64) {
+ // This is based on Peter John Acklam's inverse normal CDF
+ // algorithm: http://home.online.no/~pjacklam/notes/invnorm/
+ const (
+ a1 = -3.969683028665376e+01
+ a2 = 2.209460984245205e+02
+ a3 = -2.759285104469687e+02
+ a4 = 1.383577518672690e+02
+ a5 = -3.066479806614716e+01
+ a6 = 2.506628277459239e+00
+
+ b1 = -5.447609879822406e+01
+ b2 = 1.615858368580409e+02
+ b3 = -1.556989798598866e+02
+ b4 = 6.680131188771972e+01
+ b5 = -1.328068155288572e+01
+
+ c1 = -7.784894002430293e-03
+ c2 = -3.223964580411365e-01
+ c3 = -2.400758277161838e+00
+ c4 = -2.549732539343734e+00
+ c5 = 4.374664141464968e+00
+ c6 = 2.938163982698783e+00
+
+ d1 = 7.784695709041462e-03
+ d2 = 3.224671290700398e-01
+ d3 = 2.445134137142996e+00
+ d4 = 3.754408661907416e+00
+
+ plow = 0.02425
+ phigh = 1 - plow
+ )
+
+ if p < 0 || p > 1 {
+ return nan
+ } else if p == 0 {
+ return -inf
+ } else if p == 1 {
+ return inf
+ }
+
+ if p < plow {
+ // Rational approximation for lower region.
+ q := math.Sqrt(-2 * math.Log(p))
+ x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) /
+ ((((d1*q+d2)*q+d3)*q+d4)*q + 1)
+ } else if phigh < p {
+ // Rational approximation for upper region.
+ q := math.Sqrt(-2 * math.Log(1-p))
+ x = -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q + c6) /
+ ((((d1*q+d2)*q+d3)*q+d4)*q + 1)
+ } else {
+ // Rational approximation for central region.
+ q := p - 0.5
+ r := q * q
+ x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r + a6) * q /
+ (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r + 1)
+ }
+
+ // Refine approximation.
+ e := 0.5*math.Erfc(-x/math.Sqrt2) - p
+ u := e * math.Sqrt(2*math.Pi) * math.Exp(x*x/2)
+ x = x - u/(1+x*u/2)
+
+ // Adjust from standard normal.
+ return x*n.Sigma + n.Mu
+}
+
+func (n NormalDist) Rand(r *rand.Rand) float64 {
+ var x float64
+ if r == nil {
+ x = rand.NormFloat64()
+ } else {
+ x = r.NormFloat64()
+ }
+ return x*n.Sigma + n.Mu
+}
+
+func (n NormalDist) Bounds() (float64, float64) {
+ const stddevs = 3
+ return n.Mu - stddevs*n.Sigma, n.Mu + stddevs*n.Sigma
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/package.go b/vendor/github.com/aclements/go-moremath/stats/package.go
new file mode 100644
index 0000000..644b399
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/package.go
@@ -0,0 +1,25 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package stats implements several statistical distributions,
+// hypothesis tests, and functions for descriptive statistics.
+//
+// Currently stats is fairly small, but for what it does implement, it
+// focuses on high quality, fast implementations with good, idiomatic
+// Go APIs.
+package stats // import "github.com/aclements/go-moremath/stats"
+
+import (
+ "errors"
+ "math"
+)
+
+var inf = math.Inf(1)
+var nan = math.NaN()
+
+// TODO: Put all errors in the same place and maybe unify them.
+
+var (
+ ErrSamplesEqual = errors.New("all samples are equal")
+)
diff --git a/vendor/github.com/aclements/go-moremath/stats/sample.go b/vendor/github.com/aclements/go-moremath/stats/sample.go
new file mode 100644
index 0000000..0b5d23e
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/sample.go
@@ -0,0 +1,335 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+ "sort"
+
+ "github.com/aclements/go-moremath/vec"
+)
+
+// Sample is a collection of possibly weighted data points.
+type Sample struct {
+ // Xs is the slice of sample values.
+ Xs []float64
+
+ // Weights[i] is the weight of sample Xs[i]. If Weights is
+ // nil, all Xs have weight 1. Weights must have the same
+ // length of Xs and all values must be non-negative.
+ Weights []float64
+
+ // Sorted indicates that Xs is sorted in ascending order.
+ Sorted bool
+}
+
+// Bounds returns the minimum and maximum values of xs.
+func Bounds(xs []float64) (min float64, max float64) {
+ if len(xs) == 0 {
+ return math.NaN(), math.NaN()
+ }
+ min, max = xs[0], xs[0]
+ for _, x := range xs {
+ if x < min {
+ min = x
+ }
+ if x > max {
+ max = x
+ }
+ }
+ return
+}
+
+// Bounds returns the minimum and maximum values of the Sample.
+//
+// If the Sample is weighted, this ignores samples with zero weight.
+//
+// This is constant time if s.Sorted and there are no zero-weighted
+// values.
+func (s Sample) Bounds() (min float64, max float64) {
+ if len(s.Xs) == 0 || (!s.Sorted && s.Weights == nil) {
+ return Bounds(s.Xs)
+ }
+
+ if s.Sorted {
+ if s.Weights == nil {
+ return s.Xs[0], s.Xs[len(s.Xs)-1]
+ }
+ min, max = math.NaN(), math.NaN()
+ for i, w := range s.Weights {
+ if w != 0 {
+ min = s.Xs[i]
+ break
+ }
+ }
+ if math.IsNaN(min) {
+ return
+ }
+ for i := range s.Weights {
+ if s.Weights[len(s.Weights)-i-1] != 0 {
+ max = s.Xs[len(s.Weights)-i-1]
+ break
+ }
+ }
+ } else {
+ min, max = math.Inf(1), math.Inf(-1)
+ for i, x := range s.Xs {
+ w := s.Weights[i]
+ if x < min && w != 0 {
+ min = x
+ }
+ if x > max && w != 0 {
+ max = x
+ }
+ }
+ if math.IsInf(min, 0) {
+ min, max = math.NaN(), math.NaN()
+ }
+ }
+ return
+}
+
+// Sum returns the (possibly weighted) sum of the Sample.
+func (s Sample) Sum() float64 {
+ if s.Weights == nil {
+ return vec.Sum(s.Xs)
+ }
+ sum := 0.0
+ for i, x := range s.Xs {
+ sum += x * s.Weights[i]
+ }
+ return sum
+}
+
+// Weight returns the total weight of the Sasmple.
+func (s Sample) Weight() float64 {
+ if s.Weights == nil {
+ return float64(len(s.Xs))
+ }
+ return vec.Sum(s.Weights)
+}
+
+// Mean returns the arithmetic mean of xs.
+func Mean(xs []float64) float64 {
+ if len(xs) == 0 {
+ return math.NaN()
+ }
+ m := 0.0
+ for i, x := range xs {
+ m += (x - m) / float64(i+1)
+ }
+ return m
+}
+
+// Mean returns the arithmetic mean of the Sample.
+func (s Sample) Mean() float64 {
+ if len(s.Xs) == 0 || s.Weights == nil {
+ return Mean(s.Xs)
+ }
+
+ m, wsum := 0.0, 0.0
+ for i, x := range s.Xs {
+ // Use weighted incremental mean:
+ // m_i = (1 - w_i/wsum_i) * m_(i-1) + (w_i/wsum_i) * x_i
+ // = m_(i-1) + (x_i - m_(i-1)) * (w_i/wsum_i)
+ w := s.Weights[i]
+ wsum += w
+ m += (x - m) * w / wsum
+ }
+ return m
+}
+
+// GeoMean returns the geometric mean of xs. xs must be positive.
+func GeoMean(xs []float64) float64 {
+ if len(xs) == 0 {
+ return math.NaN()
+ }
+ m := 0.0
+ for i, x := range xs {
+ if x <= 0 {
+ return math.NaN()
+ }
+ lx := math.Log(x)
+ m += (lx - m) / float64(i+1)
+ }
+ return math.Exp(m)
+}
+
+// GeoMean returns the geometric mean of the Sample. All samples
+// values must be positive.
+func (s Sample) GeoMean() float64 {
+ if len(s.Xs) == 0 || s.Weights == nil {
+ return GeoMean(s.Xs)
+ }
+
+ m, wsum := 0.0, 0.0
+ for i, x := range s.Xs {
+ w := s.Weights[i]
+ wsum += w
+ lx := math.Log(x)
+ m += (lx - m) * w / wsum
+ }
+ return math.Exp(m)
+}
+
+// Variance returns the sample variance of xs.
+func Variance(xs []float64) float64 {
+ if len(xs) == 0 {
+ return math.NaN()
+ } else if len(xs) <= 1 {
+ return 0
+ }
+
+ // Based on Wikipedia's presentation of Welford 1962
+ // (http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm).
+ // This is more numerically stable than the standard two-pass
+ // formula and not prone to massive cancellation.
+ mean, M2 := 0.0, 0.0
+ for n, x := range xs {
+ delta := x - mean
+ mean += delta / float64(n+1)
+ M2 += delta * (x - mean)
+ }
+ return M2 / float64(len(xs)-1)
+}
+
+func (s Sample) Variance() float64 {
+ if len(s.Xs) == 0 || s.Weights == nil {
+ return Variance(s.Xs)
+ }
+ // TODO(austin)
+ panic("Weighted Variance not implemented")
+}
+
+// StdDev returns the sample standard deviation of xs.
+func StdDev(xs []float64) float64 {
+ return math.Sqrt(Variance(xs))
+}
+
+// StdDev returns the sample standard deviation of the Sample.
+func (s Sample) StdDev() float64 {
+ if len(s.Xs) == 0 || s.Weights == nil {
+ return StdDev(s.Xs)
+ }
+ // TODO(austin)
+ panic("Weighted StdDev not implemented")
+}
+
+// Quantile returns the sample value X at which q*weight of the sample
+// is <= X. This uses interpolation method R8 from Hyndman and Fan
+// (1996).
+//
+// q will be capped to the range [0, 1]. If len(xs) == 0 or all
+// weights are 0, returns NaN.
+//
+// Quantile(0.5) is the median. Quantile(0.25) and Quantile(0.75) are
+// the first and third quartiles, respectively. Quantile(P/100) is the
+// P'th percentile.
+//
+// This is constant time if s.Sorted and s.Weights == nil.
+func (s Sample) Quantile(q float64) float64 {
+ if len(s.Xs) == 0 {
+ return math.NaN()
+ } else if q <= 0 {
+ min, _ := s.Bounds()
+ return min
+ } else if q >= 1 {
+ _, max := s.Bounds()
+ return max
+ }
+
+ if !s.Sorted {
+ // TODO(austin) Use select algorithm instead
+ s = *s.Copy().Sort()
+ }
+
+ if s.Weights == nil {
+ N := float64(len(s.Xs))
+ //n := q * (N + 1) // R6
+ n := 1/3.0 + q*(N+1/3.0) // R8
+ kf, frac := math.Modf(n)
+ k := int(kf)
+ if k <= 0 {
+ return s.Xs[0]
+ } else if k >= len(s.Xs) {
+ return s.Xs[len(s.Xs)-1]
+ }
+ return s.Xs[k-1] + frac*(s.Xs[k]-s.Xs[k-1])
+ } else {
+ // TODO(austin): Implement interpolation
+
+ target := s.Weight() * q
+
+ // TODO(austin) If we had cumulative weights, we could
+ // do this in log time.
+ for i, weight := range s.Weights {
+ target -= weight
+ if target < 0 {
+ return s.Xs[i]
+ }
+ }
+ return s.Xs[len(s.Xs)-1]
+ }
+}
+
+// IQR returns the interquartile range of the Sample.
+//
+// This is constant time if s.Sorted and s.Weights == nil.
+func (s Sample) IQR() float64 {
+ if !s.Sorted {
+ s = *s.Copy().Sort()
+ }
+ return s.Quantile(0.75) - s.Quantile(0.25)
+}
+
+type sampleSorter struct {
+ xs []float64
+ weights []float64
+}
+
+func (p *sampleSorter) Len() int {
+ return len(p.xs)
+}
+
+func (p *sampleSorter) Less(i, j int) bool {
+ return p.xs[i] < p.xs[j]
+}
+
+func (p *sampleSorter) Swap(i, j int) {
+ p.xs[i], p.xs[j] = p.xs[j], p.xs[i]
+ p.weights[i], p.weights[j] = p.weights[j], p.weights[i]
+}
+
+// Sort sorts the samples in place in s and returns s.
+//
+// A sorted sample improves the performance of some algorithms.
+func (s *Sample) Sort() *Sample {
+ if s.Sorted || sort.Float64sAreSorted(s.Xs) {
+ // All set
+ } else if s.Weights == nil {
+ sort.Float64s(s.Xs)
+ } else {
+ sort.Sort(&sampleSorter{s.Xs, s.Weights})
+ }
+ s.Sorted = true
+ return s
+}
+
+// Copy returns a copy of the Sample.
+//
+// The returned Sample shares no data with the original, so they can
+// be modified (for example, sorted) independently.
+func (s Sample) Copy() *Sample {
+ xs := make([]float64, len(s.Xs))
+ copy(xs, s.Xs)
+
+ weights := []float64(nil)
+ if s.Weights != nil {
+ weights = make([]float64, len(s.Weights))
+ copy(weights, s.Weights)
+ }
+
+ return &Sample{xs, weights, s.Sorted}
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/stream.go b/vendor/github.com/aclements/go-moremath/stats/stream.go
new file mode 100644
index 0000000..0deb904
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/stream.go
@@ -0,0 +1,100 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "fmt"
+ "math"
+)
+
+// TODO(austin) Unify more with Sample interface
+
+// StreamStats tracks basic statistics for a stream of data in O(1)
+// space.
+//
+// StreamStats should be initialized to its zero value.
+type StreamStats struct {
+ Count uint
+ Total, Min, Max float64
+
+ // Numerically stable online mean
+ mean float64
+ meanOfSquares float64
+
+ // Online variance
+ vM2 float64
+}
+
+// Add updates s's statistics with sample value x.
+func (s *StreamStats) Add(x float64) {
+ s.Total += x
+ if s.Count == 0 {
+ s.Min, s.Max = x, x
+ } else {
+ if x < s.Min {
+ s.Min = x
+ }
+ if x > s.Max {
+ s.Max = x
+ }
+ }
+ s.Count++
+
+ // Update online mean, mean of squares, and variance. Online
+ // variance based on Wikipedia's presentation ("Algorithms for
+ // calculating variance") of Knuth's formulation of Welford
+ // 1962.
+ delta := x - s.mean
+ s.mean += delta / float64(s.Count)
+ s.meanOfSquares += (x*x - s.meanOfSquares) / float64(s.Count)
+ s.vM2 += delta * (x - s.mean)
+}
+
+func (s *StreamStats) Weight() float64 {
+ return float64(s.Count)
+}
+
+func (s *StreamStats) Mean() float64 {
+ return s.mean
+}
+
+func (s *StreamStats) Variance() float64 {
+ return s.vM2 / float64(s.Count-1)
+}
+
+func (s *StreamStats) StdDev() float64 {
+ return math.Sqrt(s.Variance())
+}
+
+func (s *StreamStats) RMS() float64 {
+ return math.Sqrt(s.meanOfSquares)
+}
+
+// Combine updates s's statistics as if all samples added to o were
+// added to s.
+func (s *StreamStats) Combine(o *StreamStats) {
+ count := s.Count + o.Count
+
+ // Compute combined online variance statistics
+ delta := o.mean - s.mean
+ mean := s.mean + delta*float64(o.Count)/float64(count)
+ vM2 := s.vM2 + o.vM2 + delta*delta*float64(s.Count)*float64(o.Count)/float64(count)
+
+ s.Count = count
+ s.Total += o.Total
+ if o.Min < s.Min {
+ s.Min = o.Min
+ }
+ if o.Max > s.Max {
+ s.Max = o.Max
+ }
+ s.mean = mean
+ s.meanOfSquares += (o.meanOfSquares - s.meanOfSquares) * float64(o.Count) / float64(count)
+ s.vM2 = vM2
+}
+
+func (s *StreamStats) String() string {
+ return fmt.Sprintf("Count=%d Total=%g Min=%g Mean=%g RMS=%g Max=%g StdDev=%g", s.Count, s.Total, s.Min, s.Mean(), s.RMS(), s.Max, s.StdDev())
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/tdist.go b/vendor/github.com/aclements/go-moremath/stats/tdist.go
new file mode 100644
index 0000000..29bbb1a
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/tdist.go
@@ -0,0 +1,42 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+
+ "github.com/aclements/go-moremath/mathx"
+)
+
+// A TDist is a Student's t-distribution with V degrees of freedom.
+type TDist struct {
+ V float64
+}
+
+func lgamma(x float64) float64 {
+ y, _ := math.Lgamma(x)
+ return y
+}
+
+func (t TDist) PDF(x float64) float64 {
+ return math.Exp(lgamma((t.V+1)/2)-lgamma(t.V/2)) /
+ math.Sqrt(t.V*math.Pi) * math.Pow(1+(x*x)/t.V, -(t.V+1)/2)
+}
+
+func (t TDist) CDF(x float64) float64 {
+ if x == 0 {
+ return 0.5
+ } else if x > 0 {
+ return 1 - 0.5*mathx.BetaInc(t.V/(t.V+x*x), t.V/2, 0.5)
+ } else if x < 0 {
+ return 1 - t.CDF(-x)
+ } else {
+ return math.NaN()
+ }
+}
+
+func (t TDist) Bounds() (float64, float64) {
+ return -4, 4
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/ttest.go b/vendor/github.com/aclements/go-moremath/stats/ttest.go
new file mode 100644
index 0000000..8742298
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/ttest.go
@@ -0,0 +1,147 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "errors"
+ "math"
+)
+
+// A TTestResult is the result of a t-test.
+type TTestResult struct {
+ // N1 and N2 are the sizes of the input samples. For a
+ // one-sample t-test, N2 is 0.
+ N1, N2 int
+
+ // T is the value of the t-statistic for this t-test.
+ T float64
+
+ // DoF is the degrees of freedom for this t-test.
+ DoF float64
+
+ // AltHypothesis specifies the alternative hypothesis tested
+ // by this test against the null hypothesis that there is no
+ // difference in the means of the samples.
+ AltHypothesis LocationHypothesis
+
+ // P is p-value for this t-test for the given null hypothesis.
+ P float64
+}
+
+func newTTestResult(n1, n2 int, t, dof float64, alt LocationHypothesis) *TTestResult {
+ dist := TDist{dof}
+ var p float64
+ switch alt {
+ case LocationDiffers:
+ p = 2 * (1 - dist.CDF(math.Abs(t)))
+ case LocationLess:
+ p = dist.CDF(t)
+ case LocationGreater:
+ p = 1 - dist.CDF(t)
+ }
+ return &TTestResult{N1: n1, N2: n2, T: t, DoF: dof, AltHypothesis: alt, P: p}
+}
+
+// A TTestSample is a sample that can be used for a one or two sample
+// t-test.
+type TTestSample interface {
+ Weight() float64
+ Mean() float64
+ Variance() float64
+}
+
+var (
+ ErrSampleSize = errors.New("sample is too small")
+ ErrZeroVariance = errors.New("sample has zero variance")
+ ErrMismatchedSamples = errors.New("samples have different lengths")
+)
+
+// TwoSampleTTest performs a two-sample (unpaired) Student's t-test on
+// samples x1 and x2. This is a test of the null hypothesis that x1
+// and x2 are drawn from populations with equal means. It assumes x1
+// and x2 are independent samples, that the distributions have equal
+// variance, and that the populations are normally distributed.
+func TwoSampleTTest(x1, x2 TTestSample, alt LocationHypothesis) (*TTestResult, error) {
+ n1, n2 := x1.Weight(), x2.Weight()
+ if n1 == 0 || n2 == 0 {
+ return nil, ErrSampleSize
+ }
+ v1, v2 := x1.Variance(), x2.Variance()
+ if v1 == 0 && v2 == 0 {
+ return nil, ErrZeroVariance
+ }
+
+ dof := n1 + n2 - 2
+ v12 := ((n1-1)*v1 + (n2-1)*v2) / dof
+ t := (x1.Mean() - x2.Mean()) / math.Sqrt(v12*(1/n1+1/n2))
+ return newTTestResult(int(n1), int(n2), t, dof, alt), nil
+}
+
+// TwoSampleWelchTTest performs a two-sample (unpaired) Welch's t-test
+// on samples x1 and x2. This is like TwoSampleTTest, but does not
+// assume the distributions have equal variance.
+func TwoSampleWelchTTest(x1, x2 TTestSample, alt LocationHypothesis) (*TTestResult, error) {
+ n1, n2 := x1.Weight(), x2.Weight()
+ if n1 <= 1 || n2 <= 1 {
+ // TODO: Can we still do this with n == 1?
+ return nil, ErrSampleSize
+ }
+ v1, v2 := x1.Variance(), x2.Variance()
+ if v1 == 0 && v2 == 0 {
+ return nil, ErrZeroVariance
+ }
+
+ dof := math.Pow(v1/n1+v2/n2, 2) /
+ (math.Pow(v1/n1, 2)/(n1-1) + math.Pow(v2/n2, 2)/(n2-1))
+ s := math.Sqrt(v1/n1 + v2/n2)
+ t := (x1.Mean() - x2.Mean()) / s
+ return newTTestResult(int(n1), int(n2), t, dof, alt), nil
+}
+
+// PairedTTest performs a two-sample paired t-test on samples x1 and
+// x2. If μ0 is non-zero, this tests if the average of the difference
+// is significantly different from μ0. If x1 and x2 are identical,
+// this returns nil.
+func PairedTTest(x1, x2 []float64, μ0 float64, alt LocationHypothesis) (*TTestResult, error) {
+ if len(x1) != len(x2) {
+ return nil, ErrMismatchedSamples
+ }
+ if len(x1) <= 1 {
+ // TODO: Can we still do this with n == 1?
+ return nil, ErrSampleSize
+ }
+
+ dof := float64(len(x1) - 1)
+
+ diff := make([]float64, len(x1))
+ for i := range x1 {
+ diff[i] = x1[i] - x2[i]
+ }
+ sd := StdDev(diff)
+ if sd == 0 {
+ // TODO: Can we still do the test?
+ return nil, ErrZeroVariance
+ }
+ t := (Mean(diff) - μ0) * math.Sqrt(float64(len(x1))) / sd
+ return newTTestResult(len(x1), len(x2), t, dof, alt), nil
+}
+
+// OneSampleTTest performs a one-sample t-test on sample x. This tests
+// the null hypothesis that the population mean is equal to μ0. This
+// assumes the distribution of the population of sample means is
+// normal.
+func OneSampleTTest(x TTestSample, μ0 float64, alt LocationHypothesis) (*TTestResult, error) {
+ n, v := x.Weight(), x.Variance()
+ if n == 0 {
+ return nil, ErrSampleSize
+ }
+ if v == 0 {
+ // TODO: Can we still do the test?
+ return nil, ErrZeroVariance
+ }
+ dof := n - 1
+ t := (x.Mean() - μ0) * math.Sqrt(n) / math.Sqrt(v)
+ return newTTestResult(int(n), 0, t, dof, alt), nil
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/udist.go b/vendor/github.com/aclements/go-moremath/stats/udist.go
new file mode 100644
index 0000000..06c34ad
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/udist.go
@@ -0,0 +1,389 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+
+ "github.com/aclements/go-moremath/mathx"
+)
+
+// A UDist is the discrete probability distribution of the
+// Mann-Whitney U statistic for a pair of samples of sizes N1 and N2.
+//
+// The details of computing this distribution with no ties can be
+// found in Mann, Henry B.; Whitney, Donald R. (1947). "On a Test of
+// Whether one of Two Random Variables is Stochastically Larger than
+// the Other". Annals of Mathematical Statistics 18 (1): 50–60.
+// Computing this distribution in the presence of ties is described in
+// Klotz, J. H. (1966). "The Wilcoxon, Ties, and the Computer".
+// Journal of the American Statistical Association 61 (315): 772-787
+// and Cheung, Ying Kuen; Klotz, Jerome H. (1997). "The Mann Whitney
+// Wilcoxon Distribution Using Linked Lists". Statistica Sinica 7:
+// 805-813 (the former paper contains details that are glossed over in
+// the latter paper but has mathematical typesetting issues, so it's
+// easiest to get the context from the former paper and the details
+// from the latter).
+type UDist struct {
+ N1, N2 int
+
+ // T is the count of the number of ties at each rank in the
+ // input distributions. T may be nil, in which case it is
+ // assumed there are no ties (which is equivalent to an M+N
+ // slice of 1s). It must be the case that Sum(T) == M+N.
+ T []int
+}
+
+// hasTies returns true if d has any tied samples.
+func (d UDist) hasTies() bool {
+ for _, t := range d.T {
+ if t > 1 {
+ return true
+ }
+ }
+ return false
+}
+
+// p returns the p_{d.N1,d.N2} function defined by Mann, Whitney 1947
+// for values of U from 0 up to and including the U argument.
+//
+// This algorithm runs in Θ(N1*N2*U) = O(N1²N2²) time and is quite
+// fast for small values of N1 and N2. However, it does not handle ties.
+func (d UDist) p(U int) []float64 {
+ // This is a dynamic programming implementation of the
+ // recursive recurrence definition given by Mann and Whitney:
+ //
+ // p_{n,m}(U) = (n * p_{n-1,m}(U-m) + m * p_{n,m-1}(U)) / (n+m)
+ // p_{n,m}(U) = 0 if U < 0
+ // p_{0,m}(U) = p{n,0}(U) = 1 / nCr(m+n, n) if U = 0
+ // = 0 if U > 0
+ //
+ // (Note that there is a typo in the original paper. The first
+ // recursive application of p should be for U-m, not U-M.)
+ //
+ // Since p{n,m} only depends on p{n-1,m} and p{n,m-1}, we only
+ // need to store one "plane" of the three dimensional space at
+ // a time.
+ //
+ // Furthermore, p_{n,m} = p_{m,n}, so we only construct values
+ // for n <= m and obtain the rest through symmetry.
+ //
+ // We organize the computed values of p as followed:
+ //
+ // n → N
+ // m *
+ // ↓ * *
+ // * * *
+ // * * * *
+ // * * * *
+ // M * * * *
+ //
+ // where each * is a slice indexed by U. The code below
+ // computes these left-to-right, top-to-bottom, so it only
+ // stores one row of this matrix at a time. Furthermore,
+ // computing an element in a given U slice only depends on the
+ // same and smaller values of U, so we can overwrite the U
+ // slice we're computing in place as long as we start with the
+ // largest value of U. Finally, even though the recurrence
+ // depends on (n,m) above the diagonal and we use symmetry to
+ // mirror those across the diagonal to (m,n), the mirrored
+ // indexes are always available in the current row, so this
+ // mirroring does not interfere with our ability to recycle
+ // state.
+
+ N, M := d.N1, d.N2
+ if N > M {
+ N, M = M, N
+ }
+
+ memo := make([][]float64, N+1)
+ for n := range memo {
+ memo[n] = make([]float64, U+1)
+ }
+
+ for m := 0; m <= M; m++ {
+ // Compute p_{0,m}. This is zero except for U=0.
+ memo[0][0] = 1
+
+ // Compute the remainder of this row.
+ nlim := N
+ if m < nlim {
+ nlim = m
+ }
+ for n := 1; n <= nlim; n++ {
+ lp := memo[n-1] // p_{n-1,m}
+ var rp []float64
+ if n <= m-1 {
+ rp = memo[n] // p_{n,m-1}
+ } else {
+ rp = memo[m-1] // p{m-1,n} and m==n
+ }
+
+ // For a given n,m, U is at most n*m.
+ //
+ // TODO: Actually, it's at most ⌈n*m/2⌉, but
+ // then we need to use more complex symmetries
+ // in the inner loop below.
+ ulim := n * m
+ if U < ulim {
+ ulim = U
+ }
+
+ out := memo[n] // p_{n,m}
+ nplusm := float64(n + m)
+ for U1 := ulim; U1 >= 0; U1-- {
+ l := 0.0
+ if U1-m >= 0 {
+ l = float64(n) * lp[U1-m]
+ }
+ r := float64(m) * rp[U1]
+ out[U1] = (l + r) / nplusm
+ }
+ }
+ }
+ return memo[N]
+}
+
+type ukey struct {
+ n1 int // size of first sample
+ twoU int // 2*U statistic for this permutation
+}
+
+// This computes the cumulative counts of the Mann-Whitney U
+// distribution in the presence of ties using the computation from
+// Cheung, Ying Kuen; Klotz, Jerome H. (1997). "The Mann Whitney
+// Wilcoxon Distribution Using Linked Lists". Statistica Sinica 7:
+// 805-813, with much guidance from appendix L of Klotz, A
+// Computational Approach to Statistics.
+//
+// makeUmemo constructs a table memo[K][ukey{n1, 2*U}], where K is the
+// number of ranks (up to len(t)), n1 is the size of the first sample
+// (up to the n1 argument), and U is the U statistic (up to the
+// argument twoU/2). The value of an entry in the memo table is the
+// number of permutations of a sample of size n1 in a ranking with tie
+// vector t[:K] having a U statistic <= U.
+func makeUmemo(twoU, n1 int, t []int) []map[ukey]float64 {
+ // Another candidate for a fast implementation is van de Wiel,
+ // "The split-up algorithm: a fast symbolic method for
+ // computing p-values of distribution-free statistics". This
+ // is what's used by R's coin package. It's a comparatively
+ // recent publication, so it's presumably faster (or perhaps
+ // just more general) than previous techniques, but I can't
+ // get my hands on the paper.
+ //
+ // TODO: ~40% of this function's time is spent in mapassign on
+ // the assignment lines in the two loops and another ~20% in
+ // map access and iteration. Improving map behavior or
+ // replacing the maps altogether with some other constant-time
+ // structure could double performance.
+ //
+ // TODO: The worst case for this function is when there are
+ // few ties. Yet the best case overall is when there are *no*
+ // ties. Can we get the best of both worlds? Use the fast
+ // algorithm for the most part when there are few ties and mix
+ // in the general algorithm just where we need it? That's
+ // certainly possible for sub-problems where t[:k] has no
+ // ties, but that doesn't help if t[0] has a tie but nothing
+ // else does. Is it possible to rearrange the ranks without
+ // messing up our computation of the U statistic for
+ // sub-problems?
+
+ K := len(t)
+
+ // Compute a coefficients. The a slice is indexed by k (a[0]
+ // is unused).
+ a := make([]int, K+1)
+ a[1] = t[0]
+ for k := 2; k <= K; k++ {
+ a[k] = a[k-1] + t[k-2] + t[k-1]
+ }
+
+ // Create the memo table for the counts function, A. The A
+ // slice is indexed by k (A[0] is unused).
+ //
+ // In "The Mann Whitney Distribution Using Linked Lists", they
+ // use linked lists (*gasp*) for this, but within each K it's
+ // really just a memoization table, so it's faster to use a
+ // map. The outer structure is a slice indexed by k because we
+ // need to find all memo entries with certain values of k.
+ //
+ // TODO: The n1 and twoU values in the ukeys follow strict
+ // patterns. For each K value, the n1 values are every integer
+ // between two bounds. For each (K, n1) value, the twoU values
+ // are every integer multiple of a certain base between two
+ // bounds. It might be worth turning these into directly
+ // indexible slices.
+ A := make([]map[ukey]float64, K+1)
+ A[K] = map[ukey]float64{ukey{n1: n1, twoU: twoU}: 0}
+
+ // Compute memo table (k, n1, twoU) triples from high K values
+ // to low K values. This drives the recurrence relation
+ // downward to figure out all of the needed argument triples.
+ //
+ // TODO: Is it possible to generate this table bottom-up? If
+ // so, this could be a pure dynamic programming algorithm and
+ // we could discard the K dimension. We could at least store
+ // the inputs in a more compact representation that replaces
+ // the twoU dimension with an interval and a step size (as
+ // suggested by Cheung, Klotz, not that they make it at all
+ // clear *why* they're suggesting this).
+ tsum := sumint(t) // always ∑ t[0:k]
+ for k := K - 1; k >= 2; k-- {
+ tsum -= t[k]
+ A[k] = make(map[ukey]float64)
+
+ // Construct A[k] from A[k+1].
+ for A_kplus1 := range A[k+1] {
+ rkLow := maxint(0, A_kplus1.n1-tsum)
+ rkHigh := minint(A_kplus1.n1, t[k])
+ for rk := rkLow; rk <= rkHigh; rk++ {
+ twoU_k := A_kplus1.twoU - rk*(a[k+1]-2*A_kplus1.n1+rk)
+ n1_k := A_kplus1.n1 - rk
+ if twoUmin(n1_k, t[:k], a) <= twoU_k && twoU_k <= twoUmax(n1_k, t[:k], a) {
+ key := ukey{n1: n1_k, twoU: twoU_k}
+ A[k][key] = 0
+ }
+ }
+ }
+ }
+
+ // Fill counts in memo table from low K values to high K
+ // values. This unwinds the recurrence relation.
+
+ // Start with K==2 base case.
+ //
+ // TODO: Later computations depend on these, but these don't
+ // depend on anything (including each other), so if K==2, we
+ // can skip the memo table altogether.
+ if K < 2 {
+ panic("K < 2")
+ }
+ N_2 := t[0] + t[1]
+ for A_2i := range A[2] {
+ Asum := 0.0
+ r2Low := maxint(0, A_2i.n1-t[0])
+ r2High := (A_2i.twoU - A_2i.n1*(t[0]-A_2i.n1)) / N_2
+ for r2 := r2Low; r2 <= r2High; r2++ {
+ Asum += mathx.Choose(t[0], A_2i.n1-r2) *
+ mathx.Choose(t[1], r2)
+ }
+ A[2][A_2i] = Asum
+ }
+
+ // Derive counts for the rest of the memo table.
+ tsum = t[0] // always ∑ t[0:k-1]
+ for k := 3; k <= K; k++ {
+ tsum += t[k-2]
+
+ // Compute A[k] counts from A[k-1] counts.
+ for A_ki := range A[k] {
+ Asum := 0.0
+ rkLow := maxint(0, A_ki.n1-tsum)
+ rkHigh := minint(A_ki.n1, t[k-1])
+ for rk := rkLow; rk <= rkHigh; rk++ {
+ twoU_kminus1 := A_ki.twoU - rk*(a[k]-2*A_ki.n1+rk)
+ n1_kminus1 := A_ki.n1 - rk
+ x, ok := A[k-1][ukey{n1: n1_kminus1, twoU: twoU_kminus1}]
+ if !ok && twoUmax(n1_kminus1, t[:k-1], a) < twoU_kminus1 {
+ x = mathx.Choose(tsum, n1_kminus1)
+ }
+ Asum += x * mathx.Choose(t[k-1], rk)
+ }
+ A[k][A_ki] = Asum
+ }
+ }
+
+ return A
+}
+
+func twoUmin(n1 int, t, a []int) int {
+ K := len(t)
+ twoU := -n1 * n1
+ n1_k := n1
+ for k := 1; k <= K; k++ {
+ twoU_k := minint(n1_k, t[k-1])
+ twoU += twoU_k * a[k]
+ n1_k -= twoU_k
+ }
+ return twoU
+}
+
+func twoUmax(n1 int, t, a []int) int {
+ K := len(t)
+ twoU := -n1 * n1
+ n1_k := n1
+ for k := K; k > 0; k-- {
+ twoU_k := minint(n1_k, t[k-1])
+ twoU += twoU_k * a[k]
+ n1_k -= twoU_k
+ }
+ return twoU
+}
+
+func (d UDist) PMF(U float64) float64 {
+ if U < 0 || U >= 0.5+float64(d.N1*d.N2) {
+ return 0
+ }
+
+ if d.hasTies() {
+ // makeUmemo computes the CDF directly. Take its
+ // difference to get the PMF.
+ p1, ok1 := makeUmemo(int(2*U)-1, d.N1, d.T)[len(d.T)][ukey{d.N1, int(2*U) - 1}]
+ p2, ok2 := makeUmemo(int(2*U), d.N1, d.T)[len(d.T)][ukey{d.N1, int(2 * U)}]
+ if !ok1 || !ok2 {
+ panic("makeUmemo did not return expected memoization table")
+ }
+ return (p2 - p1) / mathx.Choose(d.N1+d.N2, d.N1)
+ }
+
+ // There are no ties. Use the fast algorithm. U must be integral.
+ Ui := int(math.Floor(U))
+ // TODO: Use symmetry to minimize U
+ return d.p(Ui)[Ui]
+}
+
+func (d UDist) CDF(U float64) float64 {
+ if U < 0 {
+ return 0
+ } else if U >= float64(d.N1*d.N2) {
+ return 1
+ }
+
+ if d.hasTies() {
+ // TODO: Minimize U?
+ p, ok := makeUmemo(int(2*U), d.N1, d.T)[len(d.T)][ukey{d.N1, int(2 * U)}]
+ if !ok {
+ panic("makeUmemo did not return expected memoization table")
+ }
+ return p / mathx.Choose(d.N1+d.N2, d.N1)
+ }
+
+ // There are no ties. Use the fast algorithm. U must be integral.
+ Ui := int(math.Floor(U))
+ // The distribution is symmetric around U = m * n / 2. Sum up
+ // whichever tail is smaller.
+ flip := Ui >= (d.N1*d.N2+1)/2
+ if flip {
+ Ui = d.N1*d.N2 - Ui - 1
+ }
+ pdfs := d.p(Ui)
+ p := 0.0
+ for _, pdf := range pdfs[:Ui+1] {
+ p += pdf
+ }
+ if flip {
+ p = 1 - p
+ }
+ return p
+}
+
+func (d UDist) Step() float64 {
+ return 0.5
+}
+
+func (d UDist) Bounds() (float64, float64) {
+ // TODO: More precise bounds when there are ties.
+ return 0, float64(d.N1 * d.N2)
+}
diff --git a/vendor/github.com/aclements/go-moremath/stats/utest.go b/vendor/github.com/aclements/go-moremath/stats/utest.go
new file mode 100644
index 0000000..c31e54c
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/stats/utest.go
@@ -0,0 +1,276 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package stats
+
+import (
+ "math"
+ "sort"
+
+ "github.com/aclements/go-moremath/mathx"
+)
+
+// A LocationHypothesis specifies the alternative hypothesis of a
+// location test such as a t-test or a Mann-Whitney U-test. The
+// default (zero) value is to test against the alternative hypothesis
+// that they differ.
+type LocationHypothesis int
+
+//go:generate stringer -type LocationHypothesis
+
+const (
+ // LocationLess specifies the alternative hypothesis that the
+ // location of the first sample is less than the second. This
+ // is a one-tailed test.
+ LocationLess LocationHypothesis = -1
+
+ // LocationDiffers specifies the alternative hypothesis that
+ // the locations of the two samples are not equal. This is a
+ // two-tailed test.
+ LocationDiffers LocationHypothesis = 0
+
+ // LocationGreater specifies the alternative hypothesis that
+ // the location of the first sample is greater than the
+ // second. This is a one-tailed test.
+ LocationGreater LocationHypothesis = 1
+)
+
+// A MannWhitneyUTestResult is the result of a Mann-Whitney U-test.
+type MannWhitneyUTestResult struct {
+ // N1 and N2 are the sizes of the input samples.
+ N1, N2 int
+
+ // U is the value of the Mann-Whitney U statistic for this
+ // test, generalized by counting ties as 0.5.
+ //
+ // Given the Cartesian product of the two samples, this is the
+ // number of pairs in which the value from the first sample is
+ // greater than the value of the second, plus 0.5 times the
+ // number of pairs where the values from the two samples are
+ // equal. Hence, U is always an integer multiple of 0.5 (it is
+ // a whole integer if there are no ties) in the range [0, N1*N2].
+ //
+ // U statistics always come in pairs, depending on which
+ // sample is "first". The mirror U for the other sample can be
+ // calculated as N1*N2 - U.
+ //
+ // There are many equivalent statistics with slightly
+ // different definitions. The Wilcoxon (1945) W statistic
+ // (generalized for ties) is U + (N1(N1+1))/2. It is also
+ // common to use 2U to eliminate the half steps and Smid
+ // (1956) uses N1*N2 - 2U to additionally center the
+ // distribution.
+ U float64
+
+ // AltHypothesis specifies the alternative hypothesis tested
+ // by this test against the null hypothesis that there is no
+ // difference in the locations of the samples.
+ AltHypothesis LocationHypothesis
+
+ // P is the p-value of the Mann-Whitney test for the given
+ // null hypothesis.
+ P float64
+}
+
+// MannWhitneyExactLimit gives the largest sample size for which the
+// exact U distribution will be used for the Mann-Whitney U-test.
+//
+// Using the exact distribution is necessary for small sample sizes
+// because the distribution is highly irregular. However, computing
+// the distribution for large sample sizes is both computationally
+// expensive and unnecessary because it quickly approaches a normal
+// approximation. Computing the distribution for two 50 value samples
+// takes a few milliseconds on a 2014 laptop.
+var MannWhitneyExactLimit = 50
+
+// MannWhitneyTiesExactLimit gives the largest sample size for which
+// the exact U distribution will be used for the Mann-Whitney U-test
+// in the presence of ties.
+//
+// Computing this distribution is more expensive than computing the
+// distribution without ties, so this is set lower. Computing this
+// distribution for two 25 value samples takes about ten milliseconds
+// on a 2014 laptop.
+var MannWhitneyTiesExactLimit = 25
+
+// MannWhitneyUTest performs a Mann-Whitney U-test [1,2] of the null
+// hypothesis that two samples come from the same population against
+// the alternative hypothesis that one sample tends to have larger or
+// smaller values than the other.
+//
+// This is similar to a t-test, but unlike the t-test, the
+// Mann-Whitney U-test is non-parametric (it does not assume a normal
+// distribution). It has very slightly lower efficiency than the
+// t-test on normal distributions.
+//
+// Computing the exact U distribution is expensive for large sample
+// sizes, so this uses a normal approximation for sample sizes larger
+// than MannWhitneyExactLimit if there are no ties or
+// MannWhitneyTiesExactLimit if there are ties. This normal
+// approximation uses both the tie correction and the continuity
+// correction.
+//
+// This can fail with ErrSampleSize if either sample is empty or
+// ErrSamplesEqual if all sample values are equal.
+//
+// This is also known as a Mann-Whitney-Wilcoxon test and is
+// equivalent to the Wilcoxon rank-sum test, though the Wilcoxon
+// rank-sum test differs in nomenclature.
+//
+// [1] Mann, Henry B.; Whitney, Donald R. (1947). "On a Test of
+// Whether one of Two Random Variables is Stochastically Larger than
+// the Other". Annals of Mathematical Statistics 18 (1): 50–60.
+//
+// [2] Klotz, J. H. (1966). "The Wilcoxon, Ties, and the Computer".
+// Journal of the American Statistical Association 61 (315): 772-787.
+func MannWhitneyUTest(x1, x2 []float64, alt LocationHypothesis) (*MannWhitneyUTestResult, error) {
+ n1, n2 := len(x1), len(x2)
+ if n1 == 0 || n2 == 0 {
+ return nil, ErrSampleSize
+ }
+
+ // Compute the U statistic and tie vector T.
+ x1 = append([]float64(nil), x1...)
+ x2 = append([]float64(nil), x2...)
+ sort.Float64s(x1)
+ sort.Float64s(x2)
+ merged, labels := labeledMerge(x1, x2)
+
+ R1 := 0.0
+ T, hasTies := []int{}, false
+ for i := 0; i < len(merged); {
+ rank1, nx1, v1 := i+1, 0, merged[i]
+ // Consume samples that tie this sample (including itself).
+ for ; i < len(merged) && merged[i] == v1; i++ {
+ if labels[i] == 1 {
+ nx1++
+ }
+ }
+ // Assign all tied samples the average rank of the
+ // samples, where merged[0] has rank 1.
+ if nx1 != 0 {
+ rank := float64(i+rank1) / 2
+ R1 += rank * float64(nx1)
+ }
+ T = append(T, i-rank1+1)
+ if i > rank1 {
+ hasTies = true
+ }
+ }
+ U1 := R1 - float64(n1*(n1+1))/2
+
+ // Compute the smaller of U1 and U2
+ U2 := float64(n1*n2) - U1
+ Usmall := math.Min(U1, U2)
+
+ var p float64
+ if !hasTies && n1 <= MannWhitneyExactLimit && n2 <= MannWhitneyExactLimit ||
+ hasTies && n1 <= MannWhitneyTiesExactLimit && n2 <= MannWhitneyTiesExactLimit {
+ // Use exact U distribution. U1 will be an integer.
+ if len(T) == 1 {
+ // All values are equal. Test is meaningless.
+ return nil, ErrSamplesEqual
+ }
+
+ dist := UDist{N1: n1, N2: n2, T: T}
+ switch alt {
+ case LocationDiffers:
+ if U1 == U2 {
+ // The distribution is symmetric about
+ // Usmall. Since the distribution is
+ // discrete, the CDF is discontinuous
+ // and if simply double CDF(Usmall),
+ // we'll double count the
+ // (non-infinitesimal) probability
+ // mass at Usmall. What we want is
+ // just the integral of the whole CDF,
+ // which is 1.
+ p = 1
+ } else {
+ p = dist.CDF(Usmall) * 2
+ }
+
+ case LocationLess:
+ p = dist.CDF(U1)
+
+ case LocationGreater:
+ p = 1 - dist.CDF(U1-1)
+ }
+ } else {
+ // Use normal approximation (with tie and continuity
+ // correction).
+ t := tieCorrection(T)
+ N := float64(n1 + n2)
+ μ_U := float64(n1*n2) / 2
+ σ_U := math.Sqrt(float64(n1*n2) * ((N + 1) - t/(N*(N-1))) / 12)
+ if σ_U == 0 {
+ return nil, ErrSamplesEqual
+ }
+ numer := U1 - μ_U
+ // Perform continuity correction.
+ switch alt {
+ case LocationDiffers:
+ numer -= mathx.Sign(numer) * 0.5
+ case LocationLess:
+ numer += 0.5
+ case LocationGreater:
+ numer -= 0.5
+ }
+ z := numer / σ_U
+ switch alt {
+ case LocationDiffers:
+ p = 2 * math.Min(StdNormal.CDF(z), 1-StdNormal.CDF(z))
+ case LocationLess:
+ p = StdNormal.CDF(z)
+ case LocationGreater:
+ p = 1 - StdNormal.CDF(z)
+ }
+ }
+
+ return &MannWhitneyUTestResult{N1: n1, N2: n2, U: U1,
+ AltHypothesis: alt, P: p}, nil
+}
+
+// labeledMerge merges sorted lists x1 and x2 into sorted list merged.
+// labels[i] is 1 or 2 depending on whether merged[i] is a value from
+// x1 or x2, respectively.
+func labeledMerge(x1, x2 []float64) (merged []float64, labels []byte) {
+ merged = make([]float64, len(x1)+len(x2))
+ labels = make([]byte, len(x1)+len(x2))
+
+ i, j, o := 0, 0, 0
+ for i < len(x1) && j < len(x2) {
+ if x1[i] < x2[j] {
+ merged[o] = x1[i]
+ labels[o] = 1
+ i++
+ } else {
+ merged[o] = x2[j]
+ labels[o] = 2
+ j++
+ }
+ o++
+ }
+ for ; i < len(x1); i++ {
+ merged[o] = x1[i]
+ labels[o] = 1
+ o++
+ }
+ for ; j < len(x2); j++ {
+ merged[o] = x2[j]
+ labels[o] = 2
+ o++
+ }
+ return
+}
+
+// tieCorrection computes the tie correction factor Σ_j (t_j³ - t_j)
+// where t_j is the number of ties in the j'th rank.
+func tieCorrection(ties []int) float64 {
+ t := 0
+ for _, tie := range ties {
+ t += tie*tie*tie - tie
+ }
+ return float64(t)
+}
diff --git a/vendor/github.com/aclements/go-moremath/vec/package.go b/vendor/github.com/aclements/go-moremath/vec/package.go
new file mode 100644
index 0000000..6bd2061
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/vec/package.go
@@ -0,0 +1,6 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package vec provides functions for float64 vectors.
+package vec // import "github.com/aclements/go-moremath/vec"
diff --git a/vendor/github.com/aclements/go-moremath/vec/vec.go b/vendor/github.com/aclements/go-moremath/vec/vec.go
new file mode 100644
index 0000000..970d441
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/vec/vec.go
@@ -0,0 +1,76 @@
+// Copyright 2015 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package vec
+
+import "math"
+
+// Vectorize returns a function g(xs) that applies f to each x in xs.
+//
+// f may be evaluated in parallel and in any order.
+func Vectorize(f func(float64) float64) func(xs []float64) []float64 {
+ return func(xs []float64) []float64 {
+ return Map(f, xs)
+ }
+}
+
+// Map returns f(x) for each x in xs.
+//
+// f may be evaluated in parallel and in any order.
+func Map(f func(float64) float64, xs []float64) []float64 {
+ // TODO(austin) Parallelize
+ res := make([]float64, len(xs))
+ for i, x := range xs {
+ res[i] = f(x)
+ }
+ return res
+}
+
+// Linspace returns num values spaced evenly between lo and hi,
+// inclusive. If num is 1, this returns an array consisting of lo.
+func Linspace(lo, hi float64, num int) []float64 {
+ res := make([]float64, num)
+ if num == 1 {
+ res[0] = lo
+ return res
+ }
+ for i := 0; i < num; i++ {
+ res[i] = lo + float64(i)*(hi-lo)/float64(num-1)
+ }
+ return res
+}
+
+// Logspace returns num values spaced evenly on a logarithmic scale
+// between base**lo and base**hi, inclusive.
+func Logspace(lo, hi float64, num int, base float64) []float64 {
+ res := Linspace(lo, hi, num)
+ for i, x := range res {
+ res[i] = math.Pow(base, x)
+ }
+ return res
+}
+
+// Sum returns the sum of xs.
+func Sum(xs []float64) float64 {
+ sum := 0.0
+ for _, x := range xs {
+ sum += x
+ }
+ return sum
+}
+
+// Concat returns the concatenation of its arguments. It does not
+// modify its inputs.
+func Concat(xss ...[]float64) []float64 {
+ total := 0
+ for _, xs := range xss {
+ total += len(xs)
+ }
+ out := make([]float64, total)
+ pos := 0
+ for _, xs := range xss {
+ pos += copy(out[pos:], xs)
+ }
+ return out
+}
diff --git a/vendor/vendor.json b/vendor/vendor.json
index fa0972c..12623ea 100644
--- a/vendor/vendor.json
+++ b/vendor/vendor.json
@@ -27,6 +27,54 @@
"revisionTime": "2017-02-02T22:37:47Z"
},
{
+ "checksumSHA1": "b95iQdItSVgUj7mm2lGPq6bAKuI=",
+ "path": "github.com/aclements/go-gg/generic",
+ "revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+ "revisionTime": "2016-09-21T14:28:07Z"
+ },
+ {
+ "checksumSHA1": "AnfMX9UQwGTIk1lOXvHC5neV/Ys=",
+ "path": "github.com/aclements/go-gg/generic/slice",
+ "revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+ "revisionTime": "2016-09-21T14:28:07Z"
+ },
+ {
+ "checksumSHA1": "dCOm7590kLFskNaN1qxnZqv3Uoc=",
+ "path": "github.com/aclements/go-gg/ggstat",
+ "revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+ "revisionTime": "2016-09-21T14:28:07Z"
+ },
+ {
+ "checksumSHA1": "iSkv+GYexi78UZka28H+2tdC5S0=",
+ "path": "github.com/aclements/go-gg/table",
+ "revision": "6dbb4e4fefb062aa608dfc5abc3279753d400000",
+ "revisionTime": "2016-09-21T14:28:07Z"
+ },
+ {
+ "checksumSHA1": "0R5tItJTud6pMwB/vVna0FdHsFg=",
+ "path": "github.com/aclements/go-moremath/fit",
+ "revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+ "revisionTime": "2016-10-14T18:41:02Z"
+ },
+ {
+ "checksumSHA1": "eKoKQST+S7X0jcu8IEa/URMOsbY=",
+ "path": "github.com/aclements/go-moremath/mathx",
+ "revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+ "revisionTime": "2016-10-14T18:41:02Z"
+ },
+ {
+ "checksumSHA1": "QJmKpM8L5l7aLi6zIRXBArf9mPk=",
+ "path": "github.com/aclements/go-moremath/stats",
+ "revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+ "revisionTime": "2016-10-14T18:41:02Z"
+ },
+ {
+ "checksumSHA1": "DZGrSYl93Bo71dvmtmvdxwi6rKM=",
+ "path": "github.com/aclements/go-moremath/vec",
+ "revision": "0ff62e0875ff52f29a54f430cefa15e4775c99b2",
+ "revisionTime": "2016-10-14T18:41:02Z"
+ },
+ {
"checksumSHA1": "q3Bc7JpLWBqhZ4M7oreGo34RSkc=",
"path": "github.com/golang/protobuf/proto",
"revision": "8ee79997227bf9b34611aee7946ae64735e6fd93",
