1 files changed, 223 insertions, 0 deletions
diff --git a/src/internal/trace/mud.go b/src/internal/trace/mud.go
new file mode 100644
index 0000000..8826306
--- /dev/null
+++ b/src/internal/trace/mud.go
@@ -0,0 +1,223 @@
+// 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 trace
+
+import (
+	"math"
+	"sort"
+)
+
+// mud is an updatable mutator utilization distribution.
+//
+// This is a continuous distribution of duration over mutator
+// utilization. For example, the integral from mutator utilization a
+// to b is the total duration during which the mutator utilization was
+// in the range [a, b].
+//
+// This distribution is *not* normalized (it is not a probability
+// distribution). This makes it easier to work with as it's being
+// updated.
+//
+// It is represented as the sum of scaled uniform distribution
+// functions and Dirac delta functions (which are treated as
+// degenerate uniform distributions).
+type mud struct {
+	sorted, unsorted []edge
+
+	// trackMass is the inverse cumulative sum to track as the
+	// distribution is updated.
+	trackMass float64
+	// trackBucket is the bucket in which trackMass falls. If the
+	// total mass of the distribution is < trackMass, this is
+	// len(hist).
+	trackBucket int
+	// trackSum is the cumulative sum of hist[:trackBucket]. Once
+	// trackSum >= trackMass, trackBucket must be recomputed.
+	trackSum float64
+
+	// hist is a hierarchical histogram of distribution mass.
+	hist [mudDegree]float64
+}
+
+const (
+	// mudDegree is the number of buckets in the MUD summary
+	// histogram.
+	mudDegree = 1024
+)
+
+type edge struct {
+	// At x, the function increases by y.
+	x, delta float64
+	// Additionally at x is a Dirac delta function with area dirac.
+	dirac float64
+}
+
+// add adds a uniform function over [l, r] scaled so the total weight
+// of the uniform is area. If l==r, this adds a Dirac delta function.
+func (d *mud) add(l, r, area float64) {
+	if area == 0 {
+		return
+	}
+
+	if r < l {
+		l, r = r, l
+	}
+
+	// Add the edges.
+	if l == r {
+		d.unsorted = append(d.unsorted, edge{l, 0, area})
+	} else {
+		delta := area / (r - l)
+		d.unsorted = append(d.unsorted, edge{l, delta, 0}, edge{r, -delta, 0})
+	}
+
+	// Update the histogram.
+	h := &d.hist
+	lbFloat, lf := math.Modf(l * mudDegree)
+	lb := int(lbFloat)
+	if lb >= mudDegree {
+		lb, lf = mudDegree-1, 1
+	}
+	if l == r {
+		h[lb] += area
+	} else {
+		rbFloat, rf := math.Modf(r * mudDegree)
+		rb := int(rbFloat)
+		if rb >= mudDegree {
+			rb, rf = mudDegree-1, 1
+		}
+		if lb == rb {
+			h[lb] += area
+		} else {
+			perBucket := area / (r - l) / mudDegree
+			h[lb] += perBucket * (1 - lf)
+			h[rb] += perBucket * rf
+			for i := lb + 1; i < rb; i++ {
+				h[i] += perBucket
+			}
+		}
+	}
+
+	// Update mass tracking.
+	if thresh := float64(d.trackBucket) / mudDegree; l < thresh {
+		if r < thresh {
+			d.trackSum += area
+		} else {
+			d.trackSum += area * (thresh - l) / (r - l)
+		}
+		if d.trackSum >= d.trackMass {
+			// The tracked mass now falls in a different
+			// bucket. Recompute the inverse cumulative sum.
+			d.setTrackMass(d.trackMass)
+		}
+	}
+}
+
+// setTrackMass sets the mass to track the inverse cumulative sum for.
+//
+// Specifically, mass is a cumulative duration, and the mutator
+// utilization bounds for this duration can be queried using
+// approxInvCumulativeSum.
+func (d *mud) setTrackMass(mass float64) {
+	d.trackMass = mass
+
+	// Find the bucket currently containing trackMass by computing
+	// the cumulative sum.
+	sum := 0.0
+	for i, val := range d.hist[:] {
+		newSum := sum + val
+		if newSum > mass {
+			// mass falls in bucket i.
+			d.trackBucket = i
+			d.trackSum = sum
+			return
+		}
+		sum = newSum
+	}
+	d.trackBucket = len(d.hist)
+	d.trackSum = sum
+}
+
+// approxInvCumulativeSum is like invCumulativeSum, but specifically
+// operates on the tracked mass and returns an upper and lower bound
+// approximation of the inverse cumulative sum.
+//
+// The true inverse cumulative sum will be in the range [lower, upper).
+func (d *mud) approxInvCumulativeSum() (float64, float64, bool) {
+	if d.trackBucket == len(d.hist) {
+		return math.NaN(), math.NaN(), false
+	}
+	return float64(d.trackBucket) / mudDegree, float64(d.trackBucket+1) / mudDegree, true
+}
+
+// invCumulativeSum returns x such that the integral of d from -∞ to x
+// is y. If the total weight of d is less than y, it returns the
+// maximum of the distribution and false.
+//
+// Specifically, y is a cumulative duration, and invCumulativeSum
+// returns the mutator utilization x such that at least y time has
+// been spent with mutator utilization <= x.
+func (d *mud) invCumulativeSum(y float64) (float64, bool) {
+	if len(d.sorted) == 0 && len(d.unsorted) == 0 {
+		return math.NaN(), false
+	}
+
+	// Sort edges.
+	edges := d.unsorted
+	sort.Slice(edges, func(i, j int) bool {
+		return edges[i].x < edges[j].x
+	})
+	// Merge with sorted edges.
+	d.unsorted = nil
+	if d.sorted == nil {
+		d.sorted = edges
+	} else {
+		oldSorted := d.sorted
+		newSorted := make([]edge, len(oldSorted)+len(edges))
+		i, j := 0, 0
+		for o := range newSorted {
+			if i >= len(oldSorted) {
+				copy(newSorted[o:], edges[j:])
+				break
+			} else if j >= len(edges) {
+				copy(newSorted[o:], oldSorted[i:])
+				break
+			} else if oldSorted[i].x < edges[j].x {
+				newSorted[o] = oldSorted[i]
+				i++
+			} else {
+				newSorted[o] = edges[j]
+				j++
+			}
+		}
+		d.sorted = newSorted
+	}
+
+	// Traverse edges in order computing a cumulative sum.
+	csum, rate, prevX := 0.0, 0.0, 0.0
+	for _, e := range d.sorted {
+		newCsum := csum + (e.x-prevX)*rate
+		if newCsum >= y {
+			// y was exceeded between the previous edge
+			// and this one.
+			if rate == 0 {
+				// Anywhere between prevX and
+				// e.x will do. We return e.x
+				// because that takes care of
+				// the y==0 case naturally.
+				return e.x, true
+			}
+			return (y-csum)/rate + prevX, true
+		}
+		newCsum += e.dirac
+		if newCsum >= y {
+			// y was exceeded by the Dirac delta at e.x.
+			return e.x, true
+		}
+		csum, prevX = newCsum, e.x
+		rate += e.delta
+	}
+	return prevX, false
+}