diff --git a/plotter/conrec.go b/plotter/conrec.go index 7af6f314..1d3e7277 100644 --- a/plotter/conrec.go +++ b/plotter/conrec.go @@ -32,7 +32,7 @@ type conrecLine func(i, j int, l line, height float64) // // For full details of the algorithm, see the paper at // http://paulbourke.net/papers/conrec/ -func conrec(g GridFunc, heights []float64, fn conrecLine) { +func conrec(g GridXYZ, heights []float64, fn conrecLine) { var ( p1, p2 point @@ -43,6 +43,31 @@ func conrec(g GridFunc, heights []float64, fn conrecLine) { im = [4]int{0, 1, 1, 0} jm = [4]int{0, 0, 1, 1} + // We differ from conrec.c in the assignment of a single value + // in cases (castab in conrec.c). The value of castab[1][1][1] is + // 3, but we set cases[1][1][1] to 0. + // + // axiom: When we have a section of the grid where all the + // Z values are equal, and equal to a contour height we would + // expect to have no internal segments to draw. + // + // This is covered by case g) in Paul Bourke's description of + // the CONREC algorithm (a triangle with three vertices the lie + // on the contour level). He says, "... case g above has no really + // satisfactory solution and fortunately will occur rarely with + // real arithmetic." and then goes on to show the following image: + // + // http://paulbourke.net/papers/conrec/conrec3.gif + // + // which shows case g) in the set where no edge is drawn, agreeing + // with our axiom above. + // + // However, in the iteration over sh at conrec.c +44, a triangle + // with all vertices on the plane is given sh = {0,0,0,0,0} and + // then when the switch at conrec.c +93 happens, castab resolves + // that to case 3 for all values of m. + // + // This is fixed by replacing castab/cases[1][1][1] with 0. cases = [3][3][3]int{ {{0, 0, 8}, {0, 2, 5}, {7, 6, 9}}, {{0, 3, 4}, {1, 0, 1}, {4, 3, 0}}, diff --git a/plotter/contour.go b/plotter/contour.go index f90faaf9..98f41fcf 100644 --- a/plotter/contour.go +++ b/plotter/contour.go @@ -2,11 +2,10 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -//go:generate ./list - package plotter import ( + "image/color" "math" "sort" @@ -17,29 +16,41 @@ import ( ) // Contour implements the Plotter interface, drawing -// a contour plot of the values in the GridFunc field. The -// order of rows is from high index to low index down. +// a contour plot of the values in the GridXYZ field. type Contour struct { - GridFunc GridFunc + GridXYZ GridXYZ // Levels describes the contour heights to plot. Levels []float64 - // LineStyle is the style of the contour lines. - draw.LineStyle + // LineStyles is the set of styles for contour + // lines. Line styles are are applied to each level + // in order, modulo the length of LineStyles. + LineStyles []draw.LineStyle // Palette is the color palette used to render - // the heat map. + // the heat map. If Palette is nil or has no + // defined color, the Contour LineStyle color + // is used. Palette palette.Palette + // Underflow and Overflow are colors used to draw + // contours outside the dynamic range defined + // by Min and Max. + Underflow color.Color + Overflow color.Color + // Min and Max define the dynamic range of the // heat map. Min, Max float64 } -// NewContour creates as new contour plotter for the given data, -// using the provided palette. -func NewContour(g GridFunc, levels []float64, p palette.Palette) *Contour { +// NewContour creates as new contour plotter for the given data, using +// the provided palette. If levels is nil, contours are generated for +// the 0.01, 0.05, 0.25, 0.5, 0.75, 0.95 and 0.99 quantiles. +// If g has Min and Max methods that return a float, those returned +// values are used to set the respective Contour fields. +func NewContour(g GridXYZ, levels []float64, p palette.Palette) *Contour { var min, max float64 type minMaxer interface { Min() float64 @@ -63,26 +74,63 @@ func NewContour(g GridFunc, levels []float64, p palette.Palette) *Contour { } } + if len(levels) == 0 { + levels = quantilesR7(g, defaultQuantiles) + } + return &Contour{ - GridFunc: g, - Levels: levels, - LineStyle: DefaultLineStyle, - Palette: p, - Min: min, - Max: max, + GridXYZ: g, + Levels: levels, + LineStyles: []draw.LineStyle{DefaultLineStyle}, + Palette: p, + Min: min, + Max: max, + } +} + +// Default quantiles for case where levels is not explicitly set. +var defaultQuantiles = []float64{0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99} + +// quantilesR7 returns the pth quantiles of the data in g according the the R-7 method. +// http://en.wikipedia.org/wiki/Quantile#Estimating_the_quantiles_of_a_population +func quantilesR7(g GridXYZ, p []float64) []float64 { + c, r := g.Dims() + data := make([]float64, 0, c*r) + for i := 0; i < c; i++ { + for j := 0; j < r; j++ { + if v := g.Z(i, j); !math.IsNaN(v) { + data = append(data, v) + } + } + } + sort.Float64s(data) + v := make([]float64, len(p)) + for j, q := range p { + if q == 1 { + v[j] = data[len(data)-1] + } + h := float64(len(data)-1) * q + i := int(h) + v[j] = data[i] + (h-math.Floor(h))*(data[i+1]-data[i]) } + return v } +// naive is a debugging constant. If true, Plot performs no contour path +// reconstruction, instead rendering each path segment individually. const naive = false // Plot implements the Plot method of the plot.Plotter interface. func (h *Contour) Plot(c draw.Canvas, plt *plot.Plot) { - pal := h.Palette.Colors() - if len(pal) == 0 { - panic("contour: empty palette") + if naive { + h.naivePlot(c, plt) + return } - c.SetLineStyle(h.LineStyle) + var pal []color.Color + if h.Palette != nil { + pal = h.Palette.Colors() + } trX, trY := plt.Transforms(&c) @@ -91,80 +139,131 @@ func (h *Contour) Plot(c draw.Canvas, plt *plot.Plot) { // The alternative naive approach is to draw each line segment as // conrec returns it. The integrated path approach allows graphical // optimisations and is necessary for contour fill shading. - if !naive { - cp := contourPaths(h.GridFunc, h.Levels, trX, trY) - ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0]) - if len(h.Levels) == 1 { - ps = 0 + cp := contourPaths(h.GridXYZ, h.Levels, trX, trY) + + // ps is a palette scaling factor to scale the palette uniformly + // across the given levels. This enables a discordance between the + // number of colours and the number of levels. Sorting is not + // necessary since contourPaths sorts the levels as a side effect. + ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0]) + if len(h.Levels) == 1 { + ps = 0 + } + + for i, z := range h.Levels { + if math.IsNaN(z) { + continue } + for _, pa := range cp[z] { + if isLoop(pa) { + pa.Close() + } - for _, z := range h.Levels { - if math.IsNaN(z) { - continue + style := h.LineStyles[i%len(h.LineStyles)] + var col color.Color + switch { + case z < h.Min: + col = h.Underflow + case z > h.Max: + col = h.Overflow + case len(pal) == 0: + col = style.Color + default: + col = pal[int((z-h.Levels[0])*ps+0.5)] // Apply palette scaling. } - for _, pa := range cp[z] { - if isLoop(pa) { - pa.Close() - } - col := pal[int((z-h.Levels[0])*ps+0.5)] - if col != nil { - c.SetColor(col) - c.Stroke(pa) - } + if col != nil && style.Width != 0 { + c.SetLineStyle(style) + c.SetColor(col) + c.Stroke(pa) } } - } else { - var pa vg.Path - sort.Float64s(h.Levels) - ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0]) - if len(h.Levels) == 1 { - ps = 0 - } + } +} - conrec(h.GridFunc, h.Levels, func(_, _ int, l line, z float64) { - if math.IsNaN(z) { - return - } +// naivePlot implements the a naive rendering approach for contours. +// It is here as a debugging mode since it simply draws line segments +// generated by conrec without further computation. +func (h *Contour) naivePlot(c draw.Canvas, plt *plot.Plot) { + var pal []color.Color + if h.Palette != nil { + pal = h.Palette.Colors() + } - pa = pa[:0] + trX, trY := plt.Transforms(&c) - x1, y1 := trX(l.p1.X), trY(l.p1.Y) - x2, y2 := trX(l.p2.X), trY(l.p2.Y) + // Sort levels prior to palette scaling since we can't depend on + // sorting as a side effect from calling contourPaths. + sort.Float64s(h.Levels) + // ps is a palette scaling factor to scale the palette uniformly + // across the given levels. This enables a discordance between the + // number of colours and the number of levels. + ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0]) + if len(h.Levels) == 1 { + ps = 0 + } - if !c.Contains(draw.Point{x1, y1}) || !c.Contains(draw.Point{x2, y2}) { - return - } + levelMap := make(map[float64]int) + for i, z := range h.Levels { + levelMap[z] = i + } - pa.Move(x1, y1) - pa.Line(x2, y2) - pa.Close() + // Draw each line segment as conrec generates it. + var pa vg.Path + conrec(h.GridXYZ, h.Levels, func(_, _ int, l line, z float64) { + if math.IsNaN(z) { + return + } - col := pal[int((z-h.Levels[0])*ps+0.5)] - if col != nil { - c.SetColor(col) - c.Stroke(pa) - } - }) - } + pa = pa[:0] + + x1, y1 := trX(l.p1.X), trY(l.p1.Y) + x2, y2 := trX(l.p2.X), trY(l.p2.Y) + + if !c.Contains(draw.Point{x1, y1}) || !c.Contains(draw.Point{x2, y2}) { + return + } + + pa.Move(x1, y1) + pa.Line(x2, y2) + pa.Close() + + style := h.LineStyles[levelMap[z]%len(h.LineStyles)] + var col color.Color + switch { + case z < h.Min: + col = h.Underflow + case z > h.Max: + col = h.Overflow + case len(pal) == 0: + col = style.Color + default: + col = pal[int((z-h.Levels[0])*ps+0.5)] // Apply palette scaling. + } + if col != nil && style.Width != 0 { + c.SetLineStyle(style) + c.SetColor(col) + c.Stroke(pa) + } + }) } // DataRange implements the DataRange method // of the plot.DataRanger interface. func (h *Contour) DataRange() (xmin, xmax, ymin, ymax float64) { - c, r := h.GridFunc.Dims() - return h.GridFunc.X(0), h.GridFunc.X(c - 1), h.GridFunc.Y(0), h.GridFunc.Y(r - 1) + c, r := h.GridXYZ.Dims() + return h.GridXYZ.X(0), h.GridXYZ.X(c - 1), h.GridXYZ.Y(0), h.GridXYZ.Y(r - 1) } // GlyphBoxes implements the GlyphBoxes method // of the plot.GlyphBoxer interface. func (h *Contour) GlyphBoxes(plt *plot.Plot) []plot.GlyphBox { - c, r := h.GridFunc.Dims() + c, r := h.GridXYZ.Dims() b := make([]plot.GlyphBox, 0, r*c) for i := 0; i < c; i++ { for j := 0; j < r; j++ { b = append(b, plot.GlyphBox{ - X: plt.X.Norm(h.GridFunc.X(i)), - Y: plt.Y.Norm(h.GridFunc.Y(j)), + X: plt.X.Norm(h.GridXYZ.X(i)), + Y: plt.Y.Norm(h.GridXYZ.Y(j)), Rect: draw.Rect{ Min: draw.Point{-2.5, -2.5}, Size: draw.Point{5, 5}, @@ -185,8 +284,9 @@ func isLoop(p vg.Path) bool { // contourPaths returns a collection of vg.Paths describing contour lines based // on the input data in m cut at the given levels. The trX and trY function // are coordinate transforms. The returned map contains slices of paths keyed -// on the value of the contour level. -func contourPaths(m GridFunc, levels []float64, trX, trY func(float64) vg.Length) map[float64][]vg.Path { +// on the value of the contour level. contouPaths sorts levels ascending as a +// side effect. +func contourPaths(m GridXYZ, levels []float64, trX, trY func(float64) vg.Length) map[float64][]vg.Path { sort.Float64s(levels) ends := make(map[float64]endMap) @@ -198,7 +298,7 @@ func contourPaths(m GridFunc, levels []float64, trX, trY func(float64) vg.Length // TODO(kortschak): Check that all non-loop paths have // both ends at boundary. If any end is not at a boundary - // it must have a partner near by. Find this partner and join + // it may have a partner near by. Find this partner and join // the two conts by merging the near by ends at the mean // location. This operation is done level by level to ensure // close contours of different heights are not joined. @@ -206,12 +306,11 @@ func contourPaths(m GridFunc, levels []float64, trX, trY func(float64) vg.Length // suspect that is is possible for a bi- or higher order // furcation so it may be that the path ends at middle node // of another path. This needs to be investigated. - // - // If no partner is found and not close to boundary, panic. // Excise loops from crossed paths. for c := range conts { - c.exciseLoops(conts) + // Always try to do quick excision in production if possible. + c.exciseLoops(conts, true) } // Build vg.Paths. @@ -229,11 +328,11 @@ type contourSet map[*contour]struct{} // endMap holds a working collection of available ends. type endMap map[point]*contour -// paths extends a conrecLine function to build a set of conts that represent paths along -// contour lines. It is used as the engine for a closure where ends and conts are closed -// around in a conrecLine function, and l and z are the line and height values provided -// by conrec. At the end of a conrec call, conts will contain a map keyed on the set of -// paths with each value being the height of the contour. +// paths extends a conrecLine function to build a set of contours that represent +// paths along contour lines. It is used as the engine for a closure where ends +// and conts are closed around in a conrecLine function, and l and z are the +// line and height values provided by conrec. At the end of a conrec call, +// conts will contain a map keyed on the set of paths. func paths(l line, z float64, ends map[float64]endMap, conts contourSet) { zEnds, ok := ends[z] if !ok { @@ -283,74 +382,57 @@ func paths(l line, z float64, ends map[float64]endMap, conts contourSet) { } } -// excise removes all the elements from e1 to e2 in the list -// and creates a new list holding these elements. e1 must -// precede e2 in the list, otherwise excise will panic. -// If e1 or e2 are not in the list, excise panics. -func (l *list) excise(e1, e2 *element) *list { - if e1.list != l || e2.list != l { - panic("contour: mismatched list element") - } - - var l2 list - - p := e1.prev - n := e2.next - l2.root.next = e1 - l2.root.prev = e2 - e1.prev = &l2.root - e2.next = &l2.root - p.next = n - n.prev = p - - for e := e1; e != &l2.root; e = e.next { - if e.list == nil { - panic("contour: e2 before e1") - } - e.list = &l2 - l2.len++ - l.len-- - } - return &l2 -} +// path is a set of points forming a path. +type path []point +// contour holds a set of point lying sequentially along a contour line +// at height z. type contour struct { - l *list z float64 + + // backward and forward must each always have at least one entry. + backward path + forward path } +// newContour returns a contour starting with the end points of l for the +// height z. func newContour(l line, z float64) *contour { - li := newList() - li.PushFront(l.p1) - li.PushBack(l.p2) - return &contour{l: li, z: z} + return &contour{z: z, forward: path{l.p1}, backward: path{l.p2}} } func (c *contour) path(trX, trY func(float64) vg.Length) vg.Path { var pa vg.Path - e := c.l.Front() - p := e.Value + p := c.front() pa.Move(trX(p.X), trY(p.Y)) - for e = e.Next(); e != nil; e = e.Next() { - p = e.Value + for i := len(c.backward) - 2; i >= 0; i-- { + p = c.backward[i] + pa.Line(trX(p.X), trY(p.Y)) + } + for _, p := range c.forward { pa.Line(trX(p.X), trY(p.Y)) } return pa } -func (c *contour) front() point { return c.l.Front().Value } -func (c *contour) back() point { return c.l.Back().Value } +// front returns the first point in the contour. +func (c *contour) front() point { return c.backward[len(c.backward)-1] } + +// back returns the last point in the contour +func (c *contour) back() point { return c.forward[len(c.forward)-1] } +// extend adds the line l to the contour, updating the ends map. It returns +// a boolean indicating whether the extension was successful. func (c *contour) extend(l line, ends endMap) (ok bool) { switch c.front() { case l.p1: - c.l.PushFront(l.p2) + c.backward = append(c.backward, l.p2) delete(ends, l.p1) ends[l.p2] = c return true case l.p2: - c.l.PushFront(l.p1) + c.backward = append(c.backward, l.p1) delete(ends, l.p2) ends[l.p1] = c return true @@ -358,12 +440,12 @@ func (c *contour) extend(l line, ends endMap) (ok bool) { switch c.back() { case l.p1: - c.l.PushBack(l.p2) + c.forward = append(c.forward, l.p2) delete(ends, l.p1) ends[l.p2] = c return true case l.p2: - c.l.PushBack(l.p1) + c.forward = append(c.forward, l.p1) delete(ends, l.p2) ends[l.p1] = c return true @@ -372,57 +454,205 @@ func (c *contour) extend(l line, ends endMap) (ok bool) { return false } -func (c *contour) dropFront() { c.l.Remove(c.l.Front()) } -func (c *contour) dropBack() { c.l.Remove(c.l.Back()) } -func (c *contour) connect(b *contour, ends endMap) bool { - f1 := c.front() - f2 := b.front() - b1 := c.back() - b2 := b.back() - switch { - case f1 == f2: - delete(ends, f2) - ends[b2] = c - b.dropFront() - for i, e := b.l.Len(), b.l.Front(); i > 0; i, e = i-1, e.Next() { - c.l.PushFront(e.Value) - } +// reverse reverses the order of the point in a path and returns it. +func (p path) reverse() path { + for i, j := 0, len(p)-1; i < j; i, j = i+1, j-1 { + p[i], p[j] = p[j], p[i] + } + return p +} + +// connect connects the contour b with the receiver, updating the ends map. +// It returns a boolean indicating whether the connection was successful. +func (c *contour) connect(b *contour, ends endMap) (ok bool) { + switch c.front() { + case b.front(): + delete(ends, c.front()) + ends[b.back()] = c + c.backward = append(c.backward, b.backward.reverse()[1:]...) + c.backward = append(c.backward, b.forward...) return true - case f1 == b2: - delete(ends, b2) - ends[f2] = c - b.dropBack() - c.l.PushFrontList(b.l) - ends[c.l.Front().Value] = c + case b.back(): + delete(ends, c.front()) + ends[b.front()] = c + c.backward = append(c.backward, b.forward.reverse()[1:]...) + c.backward = append(c.backward, b.backward...) return true - case b1 == f2: - delete(ends, f2) - ends[b2] = c - b.dropFront() - c.l.PushBackList(b.l) + } + + switch c.back() { + case b.front(): + delete(ends, c.back()) + ends[b.back()] = c + c.forward = append(c.forward, b.backward.reverse()[1:]...) + c.forward = append(c.forward, b.forward...) return true - case b1 == b2: - delete(ends, b2) - ends[f2] = c - b.dropBack() - for i, e := b.l.Len(), b.l.Back(); i > 0; i, e = i-1, e.Prev() { - c.l.PushBack(e.Value) - } + case b.back(): + delete(ends, c.back()) + ends[b.front()] = c + c.forward = append(c.forward, b.forward.reverse()[1:]...) + c.forward = append(c.forward, b.backward...) return true - default: - return false } + + return false +} + +// exciseLoops finds loops within the contour that do not include the +// start and end. Loops are removed from the contour and added to the +// contour set. Loop detection is performed by Johnson's algorithm for +// finding elementary cycles. +func (c *contour) exciseLoops(conts contourSet, quick bool) { + if quick { + // Find cases we can guarantee don't need + // a complete analysis. + seen := make(map[point]struct{}) + var crossOvers int + for _, p := range c.backward { + if _, ok := seen[p]; ok { + crossOvers++ + } + seen[p] = struct{}{} + } + for _, p := range c.forward[:len(c.forward)-1] { + if _, ok := seen[p]; ok { + crossOvers++ + } + seen[p] = struct{}{} + } + switch crossOvers { + case 0: + return + case 1: + c.exciseQuick(conts) + return + } + } + + wp := append(c.backward.reverse(), c.forward...) + g := graphFrom(wp) + cycles := cyclesIn(g) + if len(cycles) == 0 { + // No further work to do but clean up after ourselves. + // We should not have reached here. + c.backward.reverse() + return + } + delete(conts, c) + + // Put loops into the contour set. + for _, cyc := range cycles { + loop := wp.subpath(cyc) + conts[&contour{ + z: c.z, + backward: loop[:1:1], + forward: loop[1:], + }] = struct{}{} + } + + // Find non-loop paths and keep them. + g.remove(cycles) + paths := wp.linearPathsIn(g) + for _, p := range paths { + conts[&contour{ + z: c.z, + backward: p[:1:1], + forward: p[1:], + }] = struct{}{} + } +} + +// graphFrom returns a graph representing the point path p. +func graphFrom(p path) graph { + g := make([]set, len(p)) + seen := make(map[point]int) + for i, v := range p { + if _, ok := seen[v]; !ok { + seen[v] = i + } + } + + for i, v := range p { + e, ok := seen[v] + if ok && g[e] == nil { + g[e] = make(set) + } + if i < len(p)-1 { + g[e][seen[p[i+1]]] = struct{}{} + } + } + + return g +} + +// subpath returns a subpath given the slice of point indices +// into the path. +func (p path) subpath(i []int) path { + pa := make(path, 0, len(i)) + for _, n := range i { + pa = append(pa, p[n]) + } + return pa +} + +// linearPathsIn returns the linear paths in g created from p. +// If g contains any cycles linearPaths will panic. +func (p path) linearPathsIn(g graph) []path { + var pa []path + + var u int + for u < len(g) { + for ; u < len(g) && len(g[u]) == 0; u++ { + } + if u == len(g) { + return pa + } + var curr path + for { + if len(g[u]) == 0 { + curr = append(curr, p[u]) + pa = append(pa, curr) + if u == len(g)-1 { + return pa + } + break + } + if len(g[u]) > 1 { + panic("contour: not a linear path") + } + for v := range g[u] { + curr = append(curr, p[u]) + u = v + break + } + } + } + + return pa } -func (c *contour) exciseLoops(conts contourSet) { - f := c.l.Front() - seen := make(map[point]*element) - for e := f; e != nil; e = e.Next() { - if p, ok := seen[e.Value]; ok && e.Value != f.Value { - nl := c.l.excise(p.Next(), e) - nl.PushFront(e.Value) - conts[&contour{l: nl, z: c.z}] = struct{}{} +// exciseQuick is a heuristic approach to loop excision. It does not +// correctly identify loops in all cases, but those cases are likely +// to be rare. +func (c *contour) exciseQuick(conts contourSet) { + wp := append(c.backward.reverse(), c.forward...) + seen := make(map[point]int) + for j := 0; j < len(wp); { + p := wp[j] + if i, ok := seen[p]; ok && p != wp[0] && p != wp[len(wp)-1] { + conts[&contour{ + z: c.z, + backward: path{wp[i]}, + forward: append(path(nil), wp[i+1:j+1]...), + }] = struct{}{} + wp = append(wp[:i], wp[j:]...) + j = i + 1 + } else { + seen[p] = j + j++ } - seen[e.Value] = e } + c.backward = c.backward[:1] + c.backward[0] = wp[0] + c.forward = wp[1:] } diff --git a/plotter/contour_test.go b/plotter/contour_test.go index dacfa8c2..33a95ffd 100644 --- a/plotter/contour_test.go +++ b/plotter/contour_test.go @@ -5,6 +5,8 @@ package plotter import ( + "flag" + "fmt" "math" "math/rand" "reflect" @@ -17,6 +19,8 @@ import ( "github.com/gonum/plot/vg" ) +var visualDebug = flag.Bool("visual", false, "output images for benchmarks and test data") + type unitGrid struct{ mat64.Matrix } func (g unitGrid) Dims() (c, r int) { r, c = g.Matrix.Dims(); return c, r } @@ -36,30 +40,66 @@ func (g unitGrid) Y(r int) float64 { return float64(r) } -func TestComplexContours(t *testing.T) { - data := make([]float64, 6400) - for i := range data { - r := float64(i/80) - 40 - c := float64(i%80) - 40 - - data[i] = rand.NormFloat64()*2 + math.Hypot(r, c) +func TestHeatMapWithContour(t *testing.T) { + if !*visualDebug { + return } + m := unitGrid{mat64.NewDense(3, 4, []float64{ + 2, 1, 4, 3, + 6, 7, 2, 5, + 9, 10, 11, 12, + })} + h := NewHeatMap(m, palette.Heat(12, 1)) - m := unitGrid{mat64.NewDense(80, 80, data)} - - levels := []float64{-1, 3, 7, 9, 13, 15, 19, 23, 27, 31} + levels := []float64{1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5} c := NewContour(m, levels, palette.Rainbow(10, palette.Blue, palette.Red, 1, 1, 1)) + c.LineStyles[0].Width *= 5 plt, _ := plot.New() + + plt.Add(h) plt.Add(c) + plt.Add(NewGlyphBoxes()) plt.X.Padding = 0 plt.Y.Padding = 0 - plt.X.Max = 79.5 - plt.Y.Max = 79.5 - plt.Save(7, 7, "complex_contour.svg") + plt.X.Max = 3.5 + plt.Y.Max = 2.5 + plt.Save(7, 7, "heat.svg") } +func TestComplexContours(t *testing.T) { + if !*visualDebug { + return + } + for _, n := range []float64{0, 1, 2, 4, 8, 16, 32} { + rand.Seed(0) + data := make([]float64, 6400) + for i := range data { + r := float64(i/80) - 40 + c := float64(i%80) - 40 + + data[i] = rand.NormFloat64()*n + math.Hypot(r, c) + } + + m := unitGrid{mat64.NewDense(80, 80, data)} + + levels := []float64{-1, 3, 7, 9, 13, 15, 19, 23, 27, 31} + c := NewContour(m, levels, palette.Rainbow(10, palette.Blue, palette.Red, 1, 1, 1)) + + plt, _ := plot.New() + plt.Add(c) + + plt.X.Padding = 0 + plt.Y.Padding = 0 + plt.X.Max = 79.5 + plt.Y.Max = 79.5 + plt.Save(7, 7, fmt.Sprintf("complex_contour-%v.svg", n)) + } +} + +func unity(f float64) vg.Length { return vg.Length(f) } + func BenchmarkComplexContour0(b *testing.B) { complexContourBench(0, b) } func BenchmarkComplexContour1(b *testing.B) { complexContourBench(1, b) } func BenchmarkComplexContour2(b *testing.B) { complexContourBench(2, b) } @@ -83,8 +123,6 @@ func complexContourBench(noise float64, b *testing.B) { levels := []float64{-1, 3, 7, 9, 13, 15, 19, 23, 27, 31} - unity := func(f float64) vg.Length { return vg.Length(f) } - var p map[float64][]vg.Path b.ResetTimer() @@ -95,31 +133,6 @@ func complexContourBench(noise float64, b *testing.B) { cp = p } -func TestHeatMapWithContour(t *testing.T) { - m := unitGrid{mat64.NewDense(3, 4, []float64{ - 2, 1, 4, 3, - 6, 7, 2, 5, - 9, 10, 11, 12, - })} - h := NewHeatMap(m, palette.Heat(12, 1)) - - levels := []float64{1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5} - c := NewContour(m, levels, palette.Rainbow(10, palette.Blue, palette.Red, 1, 1, 1)) - c.Width *= 5 - - plt, _ := plot.New() - - plt.Add(h) - plt.Add(c) - plt.Add(NewGlyphBoxes()) - - plt.X.Padding = 0 - plt.Y.Padding = 0 - plt.X.Max = 3.5 - plt.Y.Max = 2.5 - plt.Save(7, 7, "heat.svg") -} - func TestContourPaths(t *testing.T) { m := unitGrid{mat64.NewDense(3, 4, []float64{ 2, 1, 4, 3, @@ -129,8 +142,6 @@ func TestContourPaths(t *testing.T) { levels := []float64{1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5} - unity := func(x float64) vg.Length { return vg.Length(x) } - var ( wantClosed = 2 gotClosed int @@ -372,53 +383,89 @@ var wantContours = map[float64][]vg.Path{ }, } -func TestExcise(t *testing.T) { - const elements = 10 - - for i := 0; i < elements; i++ { - for j := 0; j < elements; j++ { - var ( - l1 list - l2 *list - es []*element - ) - for p := 0.0; p < elements; p++ { - es = append(es, l1.PushBack(point{p, p})) - } - func() { - defer func() { - panicked := recover() != nil - if panicked != (j < i) { - t.Errorf("unexpected panic status for %d-%d: got:%t want:%t", i, j, panicked, j < i) - } - }() - l2 = l1.excise(es[i], es[j]) - }() - if i <= j { - if l1.Len() != elements-(j-i+1) { - t.Errorf("unexpected l1 length: got:%d want:%d", l1.Len(), elements-(j-i+1)) - } - if l2.Len() != j-i+1 { - t.Errorf("unexpected l2 length: got:%d want:%d", l2.Len(), j-i+1) - } +var loopTests = []struct { + c *contour - for k, e := 0, l1.Front(); e != nil; e = e.Next() { - if k == i { - k = j + 1 - } - if e.Value != es[k].Value { - t.Errorf("unexpected value: got:%v want:%v", e.Value, es[k].Value) - } - k++ - } + want []*contour +}{ + { + c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {7, 7}, {8, 8}, {9, 9}}}, + want: []*contour{ + &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {7, 7}, {8, 8}, {9, 9}}}, + }, + }, + { + c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {4, 4}, {7, 7}, {8, 8}, {9, 9}}}, + want: []*contour{ + &contour{backward: path{{4, 4}}, forward: path{{5, 5}, {6, 6}, {4, 4}}}, + &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {7, 7}, {8, 8}, {9, 9}}}, + }, + }, + { + c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {3, 3}, {7, 7}, {1, 1}, {9, 9}}}, + want: []*contour{ + &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {9, 9}}}, + &contour{backward: path{{3, 3}}, forward: path{{4, 4}, {5, 5}, {3, 3}}}, + &contour{backward: path{{1, 1}}, forward: path{{2, 2}, {3, 3}, {7, 7}, {1, 1}}}, + }, + }, + { + c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {2, 2}, {7, 7}, {2, 2}, {9, 9}}}, + want: []*contour{ + &contour{backward: path{{2, 2}}, forward: path{{7, 7}, {2, 2}}}, + &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {9, 9}}}, + &contour{backward: path{{2, 2}}, forward: path{{3, 3}, {4, 4}, {5, 5}, {2, 2}}}, + }, + }, + { + // This test is a known failing case for exciseQuick. + c: &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {3, 3}, {8, 8}, {9, 9}, {5, 5}, {10, 10}}}, + want: []*contour{ + &contour{backward: path{{5, 5}}, forward: path{{10, 10}}}, + &contour{backward: path{{0, 0}}, forward: path{{1, 1}, {2, 2}, {3, 3}}}, + &contour{backward: path{{3, 3}}, forward: path{{4, 4}, {5, 5}, {6, 6}, {3, 3}}}, + &contour{backward: path{{3, 3}}, forward: path{{8, 8}, {9, 9}, {5, 5}, {6, 6}, {3, 3}}}, + }, + }, +} - for k, e := i, l2.Front(); e != nil; e = e.Next() { - if e.Value != es[k].Value { - t.Errorf("unexpected value: got:%v want:%v", e.Value, es[k].Value) - } - k++ - } +func (c testContour) String() string { + var s string + for i, p := range c { + if i != 0 { + s += ", " + } + s += fmt.Sprintf("%v", append(p.backward.reverse(), p.forward...)) + p.backward.reverse() + } + return s +} + +func TestExciseLoops(t *testing.T) { + for _, quick := range []bool{true, false} { + for i, test := range loopTests { + gotSet := make(contourSet) + c := &contour{ + backward: append(path(nil), test.c.backward...), + forward: append(path(nil), test.c.forward...), + } + gotSet[c] = struct{}{} + c.exciseLoops(gotSet, quick) + var got []*contour + for c := range gotSet { + got = append(got, c) + } + sort.Sort(testContour(got)) + if !reflect.DeepEqual(got, test.want) { + t.Errorf("unexpected loop excision result for %d quick=%t:\n\tgot:%v\n\twant:%v", + i, quick, testContour(got), testContour(test.want)) } } } } + +type testContour []*contour + +func (c testContour) Len() int { return len(c) } +func (c testContour) Less(i, j int) bool { return len(c[i].forward) < len(c[j].forward) } +func (c testContour) Swap(i, j int) { c[i], c[j] = c[j], c[i] } diff --git a/plotter/heat.go b/plotter/heat.go index 32e74138..079ce222 100644 --- a/plotter/heat.go +++ b/plotter/heat.go @@ -14,30 +14,33 @@ import ( "github.com/gonum/plot/vg/draw" ) -// GridFunc describes three dimensional data where the X and Y +// GridXYZ describes three dimensional data where the X and Y // coordinates are arranged on a rectangular grid. -type GridFunc interface { +type GridXYZ interface { // Dims returns the dimensions of the grid. Dims() (c, r int) - // Z returns the value of a matrix element at (c, r). - // It will panic if c or r are out of bounds for the matrix. + // Z returns the value of a grid value at (c, r). + // It will panic if c or r are out of bounds for the grid. Z(c, r int) float64 // X returns the coordinate for the column at the index x. + // It will panic if c is out of bounds for the grid. X(c int) float64 + // Y returns the coordinate for the row at the index r. + // It will panic if r is out of bounds for the grid. Y(r int) float64 } // HeatMap implements the Plotter interface, drawing -// a heat map of the values in the GridFunc field. The -// order of rows is from high index to low index down. +// a heat map of the values in the GridXYZ field. type HeatMap struct { - GridFunc GridFunc + GridXYZ GridXYZ // Palette is the color palette used to render - // the heat map. + // the heat map. Palette must not be nil or + // return a zero length []color.Color. Palette palette.Palette // Underflow and Overflow are colors used to fill @@ -52,8 +55,10 @@ type HeatMap struct { } // NewHeatMap creates as new heat map plotter for the given data, -// using the provided palette. -func NewHeatMap(g GridFunc, p palette.Palette) *HeatMap { +// using the provided palette. If g has Min and Max methods that return +// a float, those returned values are used to set the respective HeatMap +// fields. +func NewHeatMap(g GridXYZ, p palette.Palette) *HeatMap { var min, max float64 type minMaxer interface { Min() float64 @@ -78,10 +83,10 @@ func NewHeatMap(g GridFunc, p palette.Palette) *HeatMap { } return &HeatMap{ - GridFunc: g, - Palette: p, - Min: min, - Max: max, + GridXYZ: g, + Palette: p, + Min: min, + Max: max, } } @@ -91,29 +96,30 @@ func (h *HeatMap) Plot(c draw.Canvas, plt *plot.Plot) { if len(pal) == 0 { panic("heatmap: empty palette") } + // ps scales the palette uniformly across the data range. ps := float64(len(pal)-1) / (h.Max - h.Min) trX, trY := plt.Transforms(&c) var pa vg.Path - cols, rows := h.GridFunc.Dims() + cols, rows := h.GridXYZ.Dims() for i := 0; i < cols; i++ { var right, left float64 switch i { case 0: - right = (h.GridFunc.X(i+1) - h.GridFunc.X(i)) / 2 + right = (h.GridXYZ.X(i+1) - h.GridXYZ.X(i)) / 2 left = -right case cols - 1: - right = (h.GridFunc.X(i) - h.GridFunc.X(i-1)) / 2 + right = (h.GridXYZ.X(i) - h.GridXYZ.X(i-1)) / 2 left = -right default: - right = (h.GridFunc.X(i+1) - h.GridFunc.X(i)) / 2 - left = -(h.GridFunc.X(i) - h.GridFunc.X(i-1)) / 2 + right = (h.GridXYZ.X(i+1) - h.GridXYZ.X(i)) / 2 + left = -(h.GridXYZ.X(i) - h.GridXYZ.X(i-1)) / 2 } for j := 0; j < rows; j++ { - v := h.GridFunc.Z(i, j) + v := h.GridXYZ.Z(i, j) if math.IsNaN(v) || math.IsInf(v, 0) { continue } @@ -123,18 +129,18 @@ func (h *HeatMap) Plot(c draw.Canvas, plt *plot.Plot) { var up, down float64 switch j { case 0: - up = (h.GridFunc.Y(j+1) - h.GridFunc.Y(j)) / 2 + up = (h.GridXYZ.Y(j+1) - h.GridXYZ.Y(j)) / 2 down = -up case rows - 1: - up = (h.GridFunc.Y(j) - h.GridFunc.Y(j-1)) / 2 + up = (h.GridXYZ.Y(j) - h.GridXYZ.Y(j-1)) / 2 down = -up default: - up = (h.GridFunc.Y(j+1) - h.GridFunc.Y(j)) / 2 - down = -(h.GridFunc.Y(j) - h.GridFunc.Y(j-1)) / 2 + up = (h.GridXYZ.Y(j+1) - h.GridXYZ.Y(j)) / 2 + down = -(h.GridXYZ.Y(j) - h.GridXYZ.Y(j-1)) / 2 } - x, y := trX(h.GridFunc.X(i)+left), trY(h.GridFunc.Y(j)+down) - dx, dy := trX(h.GridFunc.X(i)+right), trY(h.GridFunc.Y(j)+up) + x, y := trX(h.GridXYZ.X(i)+left), trY(h.GridXYZ.Y(j)+down) + dx, dy := trX(h.GridXYZ.X(i)+right), trY(h.GridXYZ.Y(j)+up) if !c.Contains(draw.Point{x, y}) || !c.Contains(draw.Point{dx, dy}) { continue @@ -153,7 +159,7 @@ func (h *HeatMap) Plot(c draw.Canvas, plt *plot.Plot) { case v > h.Max: col = h.Overflow default: - col = pal[int((v-h.Min)*ps+0.5)] + col = pal[int((v-h.Min)*ps+0.5)] // Apply palette scaling. } if col != nil { c.SetColor(col) @@ -166,22 +172,22 @@ func (h *HeatMap) Plot(c draw.Canvas, plt *plot.Plot) { // DataRange implements the DataRange method // of the plot.DataRanger interface. func (h *HeatMap) DataRange() (xmin, xmax, ymin, ymax float64) { - c, r := h.GridFunc.Dims() + c, r := h.GridXYZ.Dims() switch c { case 1: // Make a unit length when there is no neighbour. xmax = 0.5 xmin = -0.5 default: - xmax = (3*h.GridFunc.X(c-1) - h.GridFunc.X(c-2)) / 2 - xmin = (h.GridFunc.X(0) - h.GridFunc.X(1)) / 2 + xmax = (3*h.GridXYZ.X(c-1) - h.GridXYZ.X(c-2)) / 2 + xmin = (h.GridXYZ.X(0) - h.GridXYZ.X(1)) / 2 } switch r { case 1: // Make a unit length when there is no neighbour. ymax = 0.5 ymin = -0.5 default: - ymax = (3*h.GridFunc.Y(r-1) - h.GridFunc.Y(r-2)) / 2 - ymin = (h.GridFunc.Y(0) - h.GridFunc.Y(1)) / 2 + ymax = (3*h.GridXYZ.Y(r-1) - h.GridXYZ.Y(r-2)) / 2 + ymin = (h.GridXYZ.Y(0) - h.GridXYZ.Y(1)) / 2 } return xmin, xmax, ymin, ymax } @@ -189,13 +195,13 @@ func (h *HeatMap) DataRange() (xmin, xmax, ymin, ymax float64) { // GlyphBoxes implements the GlyphBoxes method // of the plot.GlyphBoxer interface. func (h *HeatMap) GlyphBoxes(plt *plot.Plot) []plot.GlyphBox { - c, r := h.GridFunc.Dims() + c, r := h.GridXYZ.Dims() b := make([]plot.GlyphBox, 0, r*c) for i := 0; i < c; i++ { for j := 0; j < r; j++ { b = append(b, plot.GlyphBox{ - X: plt.X.Norm(h.GridFunc.X(i)), - Y: plt.Y.Norm(h.GridFunc.Y(j)), + X: plt.X.Norm(h.GridXYZ.X(i)), + Y: plt.Y.Norm(h.GridXYZ.Y(j)), Rect: draw.Rect{ Min: draw.Point{-5, -5}, Size: draw.Point{10, 10}, diff --git a/plotter/johnson.go b/plotter/johnson.go new file mode 100644 index 00000000..1fc321fd --- /dev/null +++ b/plotter/johnson.go @@ -0,0 +1,277 @@ +// Copyright ©2015 The gonum 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 plotter + +// johnson implements Johnson's "Finding all the elementary +// circuits of a directed graph" algorithm. SIAM J. Comput. 4(1):1975. +// +// Comments in the johnson methods are kept in sync with the comments +// and labels from the paper. +type johnson struct { + adjacent graph // SCC adjacency list. + b []set // Johnson's "B-list". + blocked []bool + s int + + stack []int + + result [][]int +} + +// cyclesIn returns the set of elementary cycles in the graph g. +func cyclesIn(g graph) [][]int { + j := johnson{ + adjacent: g.clone(), + b: make([]set, len(g)), + blocked: make([]bool, len(g)), + } + + // len(j.adjacent) will be the order of g until Tarjan's analysis + // finds no SCC, at which point t.sccSubGraph returns nil and the + // loop breaks. + for j.s < len(j.adjacent)-1 { + // We use the previous SCC adjacency to reduce the work needed. + t := newTarjan(j.adjacent.subgraph(j.s)) + // A_k = adjacency structure of strong component K with least + // vertex in subgraph of G induced by {s, s+1, ... ,n}. + j.adjacent = t.sccSubGraph(2) // Only allow SCCs with >= 2 vertices. + if len(j.adjacent) == 0 { + break + } + + // s = least vertex in V_k + for _, v := range j.adjacent { + s := len(j.adjacent) + for n := range v { + if n < s { + s = n + } + } + if s < j.s { + j.s = s + } + } + for i, v := range j.adjacent { + if len(v) > 0 { + j.blocked[i] = false + j.b[i] = make(set) + } + } + + //L3: + _ = j.circuit(j.s) + j.s++ + } + + return j.result +} + +// circuit is the CIRCUIT sub-procedure in the paper. +func (j *johnson) circuit(v int) bool { + f := false + j.stack = append(j.stack, v) + j.blocked[v] = true + + //L1: + for w := range j.adjacent[v] { + if w == j.s { + // Output circuit composed of stack followed by s. + r := make([]int, len(j.stack)+1) + copy(r, j.stack) + r[len(r)-1] = j.s + j.result = append(j.result, r) + f = true + } else if !j.blocked[w] { + if j.circuit(w) { + f = true + } + } + } + + //L2: + if f { + j.unblock(v) + } else { + for w := range j.adjacent[v] { + j.b[w][v] = struct{}{} + } + } + j.stack = j.stack[:len(j.stack)-1] + + return f +} + +// unblock is the UNBLOCK sub-procedure in the paper. +func (j *johnson) unblock(u int) { + j.blocked[u] = false + for w := range j.b[u] { + delete(j.b[u], w) + if j.blocked[w] { + j.unblock(w) + } + } +} + +func min(a, b int) int { + if a < b { + return a + } + return b +} + +// tarjan implements Tarjan's strongly connected component finding +// algorithm. The implementation is from the pseudocode at +// +// http://en.wikipedia.org/wiki/Tarjan%27s_strongly_connected_components_algorithm +// +type tarjan struct { + g graph + + index int + indexTable []int + lowLink []int + onStack []bool + + stack []int + + sccs [][]int +} + +// newTarjan returns a tarjan with the sccs field filled with the +// strongly connected components of the directed graph g. +func newTarjan(g graph) *tarjan { + t := tarjan{ + g: g, + + indexTable: make([]int, len(g)), + lowLink: make([]int, len(g)), + onStack: make([]bool, len(g)), + } + for v := range t.g { + if t.indexTable[v] == 0 { + t.strongconnect(v) + } + } + return &t +} + +// strongconnect is the strongconnect function described in the +// wikipedia article. +func (t *tarjan) strongconnect(v int) { + // Set the depth index for v to the smallest unused index. + t.index++ + t.indexTable[v] = t.index + t.lowLink[v] = t.index + t.stack = append(t.stack, v) + t.onStack[v] = true + + // Consider successors of v. + for w := range t.g[v] { + if t.indexTable[w] == 0 { + // Successor w has not yet been visited; recur on it. + t.strongconnect(w) + t.lowLink[v] = min(t.lowLink[v], t.lowLink[w]) + } else if t.onStack[w] { + // Successor w is in stack s and hence in the current SCC. + t.lowLink[v] = min(t.lowLink[v], t.indexTable[w]) + } + } + + // If v is a root node, pop the stack and generate an SCC. + if t.lowLink[v] == t.indexTable[v] { + // Start a new strongly connected component. + var ( + scc []int + w int + ) + for { + w, t.stack = t.stack[len(t.stack)-1], t.stack[:len(t.stack)-1] + t.onStack[w] = false + // Add w to current strongly connected component. + scc = append(scc, w) + if w == v { + break + } + } + // Output the current strongly connected component. + t.sccs = append(t.sccs, scc) + } +} + +// sccSubGraph returns the graph of the tarjan's strongly connected +// components with each SCC containing at least min vertices. +// sccSubGraph returns nil if there is no SCC with at least min +// members. +func (t *tarjan) sccSubGraph(min int) graph { + if len(t.g) == 0 { + return nil + } + sub := make(graph, len(t.g)) + + var n int + for _, scc := range t.sccs { + if len(scc) < min { + continue + } + n++ + for _, u := range scc { + for _, v := range scc { + if _, ok := t.g[u][v]; ok { + if sub[u] == nil { + sub[u] = make(set) + } + sub[u][v] = struct{}{} + } + } + } + } + if n == 0 { + return nil + } + + return sub +} + +// set is an integer set. +type set map[int]struct{} + +// graph is an edge list representation of a graph. +type graph []set + +// remove deletes edges that make up the given paths from the graph. +func (g graph) remove(paths [][]int) { + for _, p := range paths { + for i, u := range p[:len(p)-1] { + delete(g[u], p[i+1]) + } + } +} + +// subgraph returns a subgraph of g induced by {s, s+1, ... , n}. The +// subgraph is destructively generated in g. +func (g graph) subgraph(s int) graph { + for u, e := range g[s:] { + for v := range e { + if v < s { + delete(g[u+s], v) + } + } + } + return g +} + +// clone returns a deep copy of the graph g. +func (g graph) clone() graph { + c := make(graph, len(g)) + for u, e := range g { + for v := range e { + if c[u] == nil { + c[u] = make(set) + } + c[u][v] = struct{}{} + } + } + return c +} diff --git a/plotter/johnson_test.go b/plotter/johnson_test.go new file mode 100644 index 00000000..a9003a0d --- /dev/null +++ b/plotter/johnson_test.go @@ -0,0 +1,394 @@ +// Copyright ©2015 The gonum 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 plotter + +import ( + "fmt" + "reflect" + "sort" + "testing" +) + +func linksTo(i ...int) set { + if len(i) == 0 { + return nil + } + s := make(set) + for _, v := range i { + s[v] = struct{}{} + } + return s +} + +func (s set) String() string { + a := make([]int, 0, len(s)) + for v := range s { + a = append(a, v) + } + sort.Ints(a) + return fmt.Sprint(a) +} + +var graphTests = []struct { + path path + g graph + + // Tarjan tests. + orderIsAmbiguous bool + wantSCCs [][]int + wantAdj graph + + // Johnson tests. + wantCycles [][]int + wantCyclePaths []path +}{ + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {7, 7}, {8, 8}, {9, 9}}, + + wantSCCs: [][]int{{9}, {8}, {7}, {6}, {5}, {4}, {3}, {2}, {1}, {0}}, + wantAdj: nil, + + wantCycles: nil, + wantCyclePaths: nil, + }, + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {4, 4}, {7, 7}, {8, 8}, {9, 9}}, + + wantSCCs: [][]int{ + {10}, {9}, {8}, + {4, 5, 6}, + {3}, {2}, {1}, {0}, + {7 /*second point{4, 4}*/}, + }, + wantAdj: graph{ + 4: linksTo(5), + 5: linksTo(6), + 6: linksTo(4), + 10: nil, + }, + + wantCycles: [][]int{ + {4, 5, 6, 4}, + }, + wantCyclePaths: []path{ + {{4, 4}, {5, 5}, {6, 6}, {4, 4}}, + }, + }, + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {4, 4}, {7, 7}, {2, 2}, {9, 9}}, + + wantSCCs: [][]int{ + {10}, + {2, 3, 4, 5, 6, 8}, + {1}, {0}, + {7 /*second point{4, 4}*/}, + {9 /*second point{2, 2}*/}, + }, + wantAdj: graph{ + 2: linksTo(3), + 3: linksTo(4), + 4: linksTo(5, 8), + 5: linksTo(6), + 6: linksTo(4), + 8: linksTo(2), + 10: nil, + }, + + wantCycles: [][]int{ + {4, 5, 6, 4}, + {2, 3, 4, 8, 2}, + }, + wantCyclePaths: []path{ + {{4, 4}, {5, 5}, {6, 6}, {4, 4}}, + {{2, 2}, {3, 3}, {4, 4}, {7, 7}, {2, 2}}, + }, + }, + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {2, 2}, {7, 7}, {2, 2}, {9, 9}}, + + wantSCCs: [][]int{{9}, + {2, 3, 4, 5, 7}, + {1}, {0}, + {6 /*second point{2, 2}*/}, + {8 /*third point{2, 2}*/}, + }, + wantAdj: graph{ + 2: linksTo(3, 7), + 3: linksTo(4), + 4: linksTo(5), + 5: linksTo(2), + 7: linksTo(2), + 9: nil, + }, + + wantCycles: [][]int{ + {2, 7, 2}, + {2, 3, 4, 5, 2}, + }, + wantCyclePaths: []path{ + {{2, 2}, {7, 7}, {2, 2}}, + {{2, 2}, {3, 3}, {4, 4}, {5, 5}, {2, 2}}, + }, + }, + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {4, 4}, {5, 5}, {6, 6}, {3, 3}, {8, 8}, {9, 9}, {5, 5}, {10, 10}}, + + wantSCCs: [][]int{ + {11}, + {3, 4, 5, 6, 8, 9}, + {2}, {1}, {0}, + {7 /*second point{4, 4}*/}, + {10 /*second point{2, 2}*/}, + }, + wantAdj: graph{ + 3: linksTo(4, 8), + 4: linksTo(5), + 5: linksTo(6), + 6: linksTo(3), + 8: linksTo(9), + 9: linksTo(5), + 11: nil, + }, + + wantCycles: [][]int{ + {3, 4, 5, 6, 3}, + {3, 8, 9, 5, 6, 3}, + }, + wantCyclePaths: []path{ + {{3, 3}, {4, 4}, {5, 5}, {6, 6}, {3, 3}}, + {{3, 3}, {8, 8}, {9, 9}, {5, 5}, {6, 6}, {3, 3}}, + }, + }, + { + path: path{{0, 0}, {1, 1}, {2, 2}, {3, 3}, {0, 0}, {5, 5}, {6, 6}, {7, 7}, {5, 5}, {9, 9}, {10, 10}, {9, 9}}, + + wantSCCs: [][]int{ + {9, 10}, + {5, 6, 7}, + {0, 1, 2, 3}, + {4 /*second point{0, 0}*/}, + {8 /*second point{5, 5}*/}, + {11 /*second point{9, 9}*/}, + }, + wantAdj: graph{ + 0: linksTo(1), + 1: linksTo(2), + 2: linksTo(3), + 3: linksTo(0), + 5: linksTo(6), + 6: linksTo(7), + 7: linksTo(5), + 9: linksTo(10), + 10: linksTo(9), + 11: nil, + }, + + wantCycles: [][]int{ + {9, 10, 9}, + {5, 6, 7, 5}, + {0, 1, 2, 3, 0}, + }, + wantCyclePaths: []path{ + {{9, 9}, {10, 10}, {9, 9}}, + {{5, 5}, {6, 6}, {7, 7}, {5, 5}}, + {{0, 0}, {1, 1}, {2, 2}, {3, 3}, {0, 0}}, + }, + }, + + { + g: graph{ + 0: linksTo(1), + 1: linksTo(2, 7), + 2: linksTo(3, 6), + 3: linksTo(4), + 4: linksTo(2, 5), + 6: linksTo(3, 5), + 7: linksTo(0, 6), + }, + + wantSCCs: [][]int{ + {5}, + {2, 3, 4, 6}, + {0, 1, 7}, + }, + wantAdj: graph{ + 0: linksTo(1), + 1: linksTo(7), + 2: linksTo(3, 6), + 3: linksTo(4), + 4: linksTo(2), + 6: linksTo(3), + 7: linksTo(0), + }, + + wantCycles: [][]int{ + {0, 1, 7, 0}, + {2, 3, 4, 2}, + {2, 6, 3, 4, 2}, + }, + }, + { + g: graph{ + 0: linksTo(1, 2, 3), + 1: linksTo(2), + 2: linksTo(3), + 3: linksTo(1), + }, + + wantSCCs: [][]int{ + {1, 2, 3}, + {0}, + }, + wantAdj: graph{ + 1: linksTo(2), + 2: linksTo(3), + 3: linksTo(1), + }, + + wantCycles: [][]int{ + {1, 2, 3, 1}, + }, + }, + { + g: graph{ + 0: linksTo(1), + 1: linksTo(0, 2), + 2: linksTo(1), + }, + + wantSCCs: [][]int{ + {0, 1, 2}, + }, + wantAdj: graph{ + 0: linksTo(1), + 1: linksTo(0, 2), + 2: linksTo(1), + }, + + wantCycles: [][]int{ + {0, 1, 0}, + {1, 2, 1}, + }, + }, + { + g: graph{ + 0: linksTo(1), + 1: linksTo(2, 3), + 2: linksTo(4, 5), + 3: linksTo(4, 5), + 4: linksTo(6), + 5: nil, + 6: nil, + }, + + orderIsAmbiguous: true, + wantSCCs: [][]int{ + // Node pairs (2, 3) and (4, 5) are not + // relatively orderable within each pair. + {6}, {5}, {4}, {3}, {2}, {1}, {0}, + }, + wantAdj: nil, + + wantCycles: nil, + }, + { + g: graph{ + 0: linksTo(1), + 1: linksTo(2, 3, 4), + 2: linksTo(0, 3), + 3: linksTo(4), + 4: linksTo(3), + }, + + orderIsAmbiguous: true, + wantSCCs: [][]int{ + // SCCs are not relatively ordable. + {3, 4}, {0, 1, 2}, + }, + wantAdj: graph{ + 0: linksTo(1), + 1: linksTo(2), + 2: linksTo(0), + 3: linksTo(4), + 4: linksTo(3), + }, + + wantCycles: [][]int{ + {3, 4, 3}, + {0, 1, 2, 0}, + }, + }, +} + +func TestTarjan(t *testing.T) { + for i, test := range graphTests { + var g graph + if test.path != nil { + g = graphFrom(test.path) + } else { + g = test.g + } + tar := newTarjan(g) + gotSCCs := tar.sccs + if test.orderIsAmbiguous { + // We lose topological order here, but that + // is not important for this use case. + sort.Sort(byComponentLengthOrStart(test.wantSCCs)) + sort.Sort(byComponentLengthOrStart(gotSCCs)) + } + // tarjan.strongconnect does range iteration over maps, + // so sort SCC members to ensure consistent ordering. + for _, scc := range gotSCCs { + sort.Ints(scc) + } + if !reflect.DeepEqual(gotSCCs, test.wantSCCs) { + t.Errorf("unexpected tarjan scc result for %d:\n\tgot:%v\n\twant:%v", i, gotSCCs, test.wantSCCs) + } + gotAdj := tar.sccSubGraph(2) + if !reflect.DeepEqual(gotAdj, test.wantAdj) { + t.Errorf("unexpected tarjan sccSubGraph(2) result for %d:\n\tgot:%#v\n\twant:%#v", i, gotAdj, test.wantAdj) + } + } +} + +func TestJohnson(t *testing.T) { + for i, test := range graphTests { + var g graph + if test.path != nil { + g = graphFrom(test.path) + } else { + g = test.g + } + gotCycles := cyclesIn(g) + // johnson.circuit does range iteration over maps, + // so sort to ensure consistent ordering. + sort.Sort(byComponentLengthOrStart(gotCycles)) + if !reflect.DeepEqual(gotCycles, test.wantCycles) { + t.Errorf("unexpected johnson result for %d:\n\tgot:%#v\n\twant:%#v", i, gotCycles, test.wantCycles) + } + + // Don't do path reconstruction tests without a path definition. + if test.path == nil { + continue + } + + // Test we reconstruct paths correctly from cycles. + var gotPaths []path + for _, pi := range gotCycles { + gotPaths = append(gotPaths, test.path.subpath(pi)) + } + if !reflect.DeepEqual(gotPaths, test.wantCyclePaths) { + t.Errorf("unexpected johnson path result for %d:\n\tgot:%#v\n\twant:%#v", i, gotPaths, test.wantCyclePaths) + } + } +} + +type byComponentLengthOrStart [][]int + +func (c byComponentLengthOrStart) Len() int { return len(c) } +func (c byComponentLengthOrStart) Less(i, j int) bool { + return len(c[i]) < len(c[j]) || (len(c[i]) == len(c[j]) && c[i][0] < c[j][0]) +} +func (c byComponentLengthOrStart) Swap(i, j int) { c[i], c[j] = c[j], c[i] } diff --git a/plotter/list b/plotter/list deleted file mode 100755 index 9d9a212e..00000000 --- a/plotter/list +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/env bash - -# Copyright ©2015 The gonum Authors. All rights reserved. -# Use of this source code is governed by a BSD-style -# license that can be found in the LICENSE file. - -cat $(go env GOROOT)/src/container/list/list.go \ -| gofmt -r 'New -> newList' \ -| gofmt -r 'List -> list' \ -| gofmt -r 'Element -> element' \ -| gofmt -r 'interface{} -> point' \ -| sed 's/package list/\npackage plotter/g' \ -> list.go diff --git a/plotter/list.go b/plotter/list.go deleted file mode 100644 index d4a076ac..00000000 --- a/plotter/list.go +++ /dev/null @@ -1,217 +0,0 @@ -// Copyright 2009 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 list implements a doubly linked list. -// -// To iterate over a list (where l is a *List): -// for e := l.Front(); e != nil; e = e.Next() { -// // do something with e.Value -// } -// - -package plotter - -// Element is an element of a linked list. -type element struct { - // Next and previous pointers in the doubly-linked list of elements. - // To simplify the implementation, internally a list l is implemented - // as a ring, such that &l.root is both the next element of the last - // list element (l.Back()) and the previous element of the first list - // element (l.Front()). - next, prev *element - - // The list to which this element belongs. - list *list - - // The value stored with this element. - Value point -} - -// Next returns the next list element or nil. -func (e *element) Next() *element { - if p := e.next; e.list != nil && p != &e.list.root { - return p - } - return nil -} - -// Prev returns the previous list element or nil. -func (e *element) Prev() *element { - if p := e.prev; e.list != nil && p != &e.list.root { - return p - } - return nil -} - -// List represents a doubly linked list. -// The zero value for List is an empty list ready to use. -type list struct { - root element // sentinel list element, only &root, root.prev, and root.next are used - len int // current list length excluding (this) sentinel element -} - -// Init initializes or clears list l. -func (l *list) Init() *list { - l.root.next = &l.root - l.root.prev = &l.root - l.len = 0 - return l -} - -// New returns an initialized list. -func newList() *list { return new(list).Init() } - -// Len returns the number of elements of list l. -// The complexity is O(1). -func (l *list) Len() int { return l.len } - -// Front returns the first element of list l or nil. -func (l *list) Front() *element { - if l.len == 0 { - return nil - } - return l.root.next -} - -// Back returns the last element of list l or nil. -func (l *list) Back() *element { - if l.len == 0 { - return nil - } - return l.root.prev -} - -// lazyInit lazily initializes a zero List value. -func (l *list) lazyInit() { - if l.root.next == nil { - l.Init() - } -} - -// insert inserts e after at, increments l.len, and returns e. -func (l *list) insert(e, at *element) *element { - n := at.next - at.next = e - e.prev = at - e.next = n - n.prev = e - e.list = l - l.len++ - return e -} - -// insertValue is a convenience wrapper for insert(&Element{Value: v}, at). -func (l *list) insertValue(v point, at *element) *element { - return l.insert(&element{Value: v}, at) -} - -// remove removes e from its list, decrements l.len, and returns e. -func (l *list) remove(e *element) *element { - e.prev.next = e.next - e.next.prev = e.prev - e.next = nil // avoid memory leaks - e.prev = nil // avoid memory leaks - e.list = nil - l.len-- - return e -} - -// Remove removes e from l if e is an element of list l. -// It returns the element value e.Value. -func (l *list) Remove(e *element) point { - if e.list == l { - // if e.list == l, l must have been initialized when e was inserted - // in l or l == nil (e is a zero Element) and l.remove will crash - l.remove(e) - } - return e.Value -} - -// PushFront inserts a new element e with value v at the front of list l and returns e. -func (l *list) PushFront(v point) *element { - l.lazyInit() - return l.insertValue(v, &l.root) -} - -// PushBack inserts a new element e with value v at the back of list l and returns e. -func (l *list) PushBack(v point) *element { - l.lazyInit() - return l.insertValue(v, l.root.prev) -} - -// InsertBefore inserts a new element e with value v immediately before mark and returns e. -// If mark is not an element of l, the list is not modified. -func (l *list) InsertBefore(v point, mark *element) *element { - if mark.list != l { - return nil - } - // see comment in List.Remove about initialization of l - return l.insertValue(v, mark.prev) -} - -// InsertAfter inserts a new element e with value v immediately after mark and returns e. -// If mark is not an element of l, the list is not modified. -func (l *list) InsertAfter(v point, mark *element) *element { - if mark.list != l { - return nil - } - // see comment in List.Remove about initialization of l - return l.insertValue(v, mark) -} - -// MoveToFront moves element e to the front of list l. -// If e is not an element of l, the list is not modified. -func (l *list) MoveToFront(e *element) { - if e.list != l || l.root.next == e { - return - } - // see comment in List.Remove about initialization of l - l.insert(l.remove(e), &l.root) -} - -// MoveToBack moves element e to the back of list l. -// If e is not an element of l, the list is not modified. -func (l *list) MoveToBack(e *element) { - if e.list != l || l.root.prev == e { - return - } - // see comment in List.Remove about initialization of l - l.insert(l.remove(e), l.root.prev) -} - -// MoveBefore moves element e to its new position before mark. -// If e or mark is not an element of l, or e == mark, the list is not modified. -func (l *list) MoveBefore(e, mark *element) { - if e.list != l || e == mark || mark.list != l { - return - } - l.insert(l.remove(e), mark.prev) -} - -// MoveAfter moves element e to its new position after mark. -// If e or mark is not an element of l, or e == mark, the list is not modified. -func (l *list) MoveAfter(e, mark *element) { - if e.list != l || e == mark || mark.list != l { - return - } - l.insert(l.remove(e), mark) -} - -// PushBackList inserts a copy of an other list at the back of list l. -// The lists l and other may be the same. -func (l *list) PushBackList(other *list) { - l.lazyInit() - for i, e := other.Len(), other.Front(); i > 0; i, e = i-1, e.Next() { - l.insertValue(e.Value, l.root.prev) - } -} - -// PushFrontList inserts a copy of an other list at the front of list l. -// The lists l and other may be the same. -func (l *list) PushFrontList(other *list) { - l.lazyInit() - for i, e := other.Len(), other.Back(); i > 0; i, e = i-1, e.Prev() { - l.insertValue(e.Value, &l.root) - } -} diff --git a/plotter/volcano b/plotter/volcano index f1220efd..a5b5065c 100755 --- a/plotter/volcano +++ b/plotter/volcano @@ -1,23 +1,27 @@ #!/usr/bin/env bash cat >volcano_example.go <> volcano_example.go '})}' -go fmt volcano_example.go \ No newline at end of file +go fmt volcano_example.go diff --git a/plotter/volcano_example.go b/plotter/volcano_example.go index 9ab173ef..85c5b87c 100644 --- a/plotter/volcano_example.go +++ b/plotter/volcano_example.go @@ -1,15 +1,181 @@ +// Generated code do not edit. Run `go generate volcano_example.go`. + // Copyright ©2015 The gonum Authors. All rights reserved. // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -// +build ignore +//go:generate ./volcano + +//+build ignore package main import ( - "fmt" + "image/color" + + "github.com/gonum/matrix/mat64" + "github.com/gonum/plot" + "github.com/gonum/plot/palette" + "github.com/gonum/plot/plotter" + "github.com/gonum/plot/vg" + "github.com/gonum/plot/vg/draw" ) +type deciGrid struct{ mat64.Matrix } + +func (g deciGrid) Dims() (c, r int) { r, c = g.Matrix.Dims(); return c, r } +func (g deciGrid) Z(c, r int) float64 { return g.Matrix.At(r, c) } +func (g deciGrid) X(c int) float64 { + _, n := g.Matrix.Dims() + if c < 0 || c >= n { + panic("index out of range") + } + return 10 * float64(c) +} +func (g deciGrid) Y(r int) float64 { + m, _ := g.Matrix.Dims() + if r < 0 || r >= m { + panic("index out of range") + } + return 10 * float64(r) +} + func main() { - fmt.Println("This file is a stub.\n\nRun `./volcano` and then `go run volcano_example.go`\n") + var levels []float64 + for l := 100.5; l < volcano.Matrix.(*mat64.Dense).Max(); l += 5 { + levels = append(levels, l) + } + c := plotter.NewContour(volcano, levels, palette.Rainbow(len(levels), (palette.Yellow+palette.Red)/2, palette.Blue, 1, 1, 1)) + quarterStyle := draw.LineStyle{ + Color: color.Black, + Width: vg.Points(0.5), + Dashes: []vg.Length{0.2, 0.4}, + } + halfStyle := draw.LineStyle{ + Color: color.Black, + Width: vg.Points(0.5), + Dashes: []vg.Length{5, 2, 1, 2}, + } + c.LineStyles = append(c.LineStyles, quarterStyle, halfStyle, quarterStyle) + + h := plotter.NewHeatMap(volcano, palette.Heat(len(levels)*2, 1)) + + p, err := plot.New() + if err != nil { + panic(err) + } + p.Title.Text = "Maunga Whau Volcano" + + p.Add(h) + p.Add(c) + + p.X.Padding = 0 + p.Y.Padding = 0 + _, p.X.Max, _, p.Y.Max = h.DataRange() + + name := "example_volcano" + + for _, ext := range []string{ + ".eps", + ".pdf", + ".svg", + ".png", + ".tiff", + ".jpg", + } { + if err := p.Save(4, 4, name+ext); err != nil { + panic(err) + } + } } + +// Data extracted from RDatasets volcano data for the Maunga Whau volcano topographic data. +var volcano = deciGrid{mat64.NewDense(87, 61, []float64{ + 100, 100, 101, 101, 101, 101, 101, 100, 100, 100, 101, 101, 102, 102, 102, 102, 103, 104, 103, 102, 101, 101, 102, 103, 104, 104, 105, 107, 107, 107, 108, 108, 110, 110, 110, 110, 110, 110, 110, 110, 108, 108, 108, 107, 107, 108, 108, 108, 108, 108, 107, 107, 107, 107, 106, 106, 105, 105, 104, 104, 103, + 101, 101, 102, 102, 102, 102, 102, 101, 101, 101, 102, 102, 103, 103, 103, 103, 104, 105, 104, 103, 102, 102, 103, 105, 106, 106, 107, 109, 110, 110, 110, 110, 111, 112, 113, 114, 116, 115, 114, 112, 110, 110, 110, 109, 108, 109, 109, 109, 109, 108, 108, 108, 108, 107, 107, 106, 106, 105, 105, 104, 104, + 102, 102, 103, 103, 103, 103, 103, 102, 102, 102, 103, 103, 104, 104, 104, 104, 105, 106, 105, 104, 104, 105, 106, 107, 108, 110, 111, 113, 114, 115, 114, 115, 116, 118, 119, 119, 121, 121, 120, 118, 116, 114, 112, 111, 110, 110, 110, 110, 109, 109, 109, 109, 108, 108, 107, 107, 106, 106, 105, 105, 104, + 103, 103, 104, 104, 104, 104, 104, 103, 103, 103, 103, 104, 104, 104, 105, 105, 106, 107, 106, 106, 106, 107, 108, 110, 111, 114, 117, 118, 117, 119, 120, 121, 122, 124, 125, 126, 127, 127, 126, 124, 122, 120, 117, 116, 113, 111, 110, 110, 110, 109, 109, 109, 109, 108, 108, 107, 107, 106, 106, 105, 105, + 104, 104, 105, 105, 105, 105, 105, 104, 104, 103, 104, 104, 105, 105, 105, 106, 107, 108, 108, 108, 109, 110, 112, 114, 115, 118, 121, 122, 121, 123, 128, 131, 129, 130, 131, 131, 132, 132, 131, 130, 128, 126, 122, 119, 115, 114, 112, 110, 110, 110, 110, 110, 109, 109, 108, 107, 107, 107, 106, 106, 105, + 105, 105, 105, 106, 106, 106, 106, 105, 105, 104, 104, 105, 105, 106, 106, 107, 109, 110, 110, 112, 113, 115, 116, 118, 119, 121, 124, 126, 126, 129, 134, 137, 137, 136, 136, 135, 136, 136, 136, 135, 133, 129, 126, 122, 118, 116, 115, 113, 111, 110, 110, 110, 110, 109, 108, 108, 108, 107, 107, 106, 106, + 105, 106, 106, 107, 107, 107, 107, 106, 106, 105, 105, 106, 106, 107, 108, 109, 111, 113, 114, 116, 118, 120, 121, 122, 123, 125, 127, 129, 130, 135, 140, 142, 142, 142, 141, 140, 140, 140, 140, 139, 137, 134, 129, 125, 121, 118, 116, 114, 112, 110, 110, 110, 111, 110, 109, 109, 108, 108, 107, 107, 106, + 106, 107, 107, 108, 108, 108, 108, 107, 107, 106, 106, 107, 108, 108, 110, 113, 115, 117, 118, 120, 122, 124, 125, 127, 128, 129, 131, 134, 135, 141, 146, 147, 146, 146, 145, 144, 144, 144, 143, 142, 141, 139, 135, 130, 126, 122, 118, 116, 114, 112, 112, 113, 112, 110, 110, 109, 109, 108, 108, 107, 106, + 107, 108, 108, 109, 109, 109, 109, 108, 108, 107, 108, 108, 110, 111, 113, 116, 118, 120, 123, 125, 127, 129, 130, 132, 134, 135, 137, 139, 142, 146, 152, 152, 151, 151, 150, 149, 148, 148, 146, 145, 143, 142, 139, 135, 131, 127, 122, 119, 117, 115, 115, 115, 114, 112, 110, 110, 109, 109, 108, 107, 107, + 108, 109, 109, 110, 110, 110, 110, 109, 109, 108, 110, 110, 113, 116, 118, 120, 122, 125, 127, 129, 133, 136, 138, 140, 141, 142, 148, 150, 151, 156, 158, 159, 158, 157, 158, 158, 154, 151, 149, 148, 146, 144, 141, 137, 134, 130, 125, 122, 120, 118, 117, 117, 115, 113, 111, 110, 110, 109, 108, 107, 107, + 109, 110, 110, 111, 111, 111, 111, 110, 110, 110, 112, 114, 118, 121, 123, 125, 127, 129, 133, 137, 141, 143, 145, 146, 148, 150, 154, 156, 159, 161, 162, 163, 164, 163, 164, 164, 160, 157, 154, 151, 149, 146, 144, 140, 137, 133, 129, 126, 124, 121, 119, 118, 116, 114, 112, 111, 110, 109, 108, 107, 106, + 110, 110, 111, 113, 112, 111, 113, 112, 112, 114, 116, 119, 121, 124, 127, 129, 133, 138, 143, 146, 149, 149, 151, 153, 154, 157, 159, 160, 163, 165, 166, 167, 168, 168, 168, 168, 166, 162, 159, 157, 154, 152, 149, 144, 140, 136, 133, 131, 128, 125, 122, 119, 117, 115, 113, 111, 110, 109, 108, 107, 106, + 110, 111, 113, 115, 114, 113, 114, 114, 115, 117, 119, 121, 124, 126, 129, 133, 140, 145, 150, 154, 155, 155, 157, 159, 161, 162, 164, 165, 167, 168, 169, 170, 172, 174, 172, 172, 171, 169, 166, 163, 161, 158, 153, 148, 143, 140, 137, 134, 131, 128, 125, 120, 118, 116, 114, 112, 110, 109, 108, 107, 105, + 111, 113, 115, 117, 116, 115, 116, 117, 117, 119, 121, 124, 126, 128, 132, 137, 143, 151, 156, 161, 161, 162, 163, 165, 166, 167, 168, 170, 171, 173, 175, 177, 179, 178, 177, 176, 176, 174, 171, 169, 165, 161, 156, 152, 148, 144, 140, 138, 135, 131, 127, 123, 119, 117, 115, 113, 111, 110, 108, 106, 105, + 114, 115, 117, 117, 117, 118, 119, 119, 120, 121, 124, 126, 128, 131, 137, 143, 150, 156, 160, 163, 165, 168, 170, 171, 172, 173, 174, 175, 177, 179, 180, 182, 183, 183, 183, 183, 180, 178, 177, 172, 168, 164, 160, 156, 152, 148, 144, 141, 138, 134, 130, 126, 121, 117, 114, 112, 110, 110, 108, 106, 104, + 116, 118, 118, 118, 120, 121, 121, 122, 122, 123, 125, 128, 130, 134, 141, 147, 152, 156, 160, 165, 168, 170, 174, 176, 179, 180, 181, 181, 182, 182, 183, 184, 186, 187, 187, 184, 184, 181, 180, 176, 172, 168, 165, 161, 157, 153, 149, 145, 142, 138, 133, 129, 125, 120, 115, 111, 110, 110, 108, 106, 104, + 118, 120, 120, 121, 122, 123, 124, 124, 125, 126, 127, 129, 132, 135, 142, 149, 153, 157, 161, 166, 170, 174, 178, 180, 182, 183, 184, 184, 185, 186, 186, 187, 189, 189, 189, 189, 189, 186, 182, 179, 175, 171, 168, 165, 162, 157, 152, 149, 145, 141, 137, 131, 125, 120, 116, 111, 110, 110, 108, 106, 104, + 120, 121, 122, 123, 124, 125, 126, 127, 127, 128, 130, 132, 134, 137, 142, 151, 155, 158, 162, 169, 172, 176, 181, 183, 184, 186, 187, 188, 189, 189, 189, 189, 190, 190, 191, 190, 190, 188, 186, 183, 180, 175, 171, 168, 165, 161, 157, 152, 149, 145, 141, 134, 127, 121, 116, 112, 110, 110, 108, 106, 104, + 120, 122, 125, 126, 126, 127, 128, 129, 130, 130, 132, 134, 136, 139, 145, 152, 157, 160, 167, 172, 175, 178, 181, 185, 186, 188, 190, 191, 192, 193, 193, 192, 192, 191, 192, 191, 191, 190, 190, 187, 184, 181, 177, 172, 169, 165, 161, 156, 152, 147, 143, 139, 131, 123, 119, 115, 111, 110, 108, 106, 105, + 121, 124, 126, 128, 129, 129, 130, 131, 132, 133, 135, 137, 139, 143, 150, 154, 159, 164, 170, 173, 176, 179, 184, 186, 189, 190, 191, 192, 193, 194, 195, 194, 193, 192, 191, 191, 191, 191, 190, 190, 188, 184, 181, 177, 173, 169, 165, 160, 155, 149, 145, 142, 136, 129, 123, 118, 114, 110, 108, 108, 107, + 122, 125, 127, 130, 130, 131, 133, 134, 135, 136, 137, 140, 143, 147, 154, 158, 162, 166, 171, 174, 177, 181, 186, 189, 190, 190, 191, 192, 191, 191, 190, 189, 188, 189, 190, 190, 191, 190, 190, 190, 189, 186, 184, 181, 177, 173, 169, 164, 158, 152, 148, 144, 140, 134, 125, 118, 115, 111, 110, 108, 107, + 122, 125, 128, 130, 132, 133, 135, 136, 137, 139, 140, 143, 147, 152, 157, 161, 164, 168, 172, 175, 179, 182, 186, 190, 190, 190, 190, 189, 187, 184, 184, 183, 182, 182, 183, 183, 183, 184, 185, 186, 187, 186, 185, 184, 181, 177, 173, 169, 163, 157, 149, 145, 141, 136, 130, 119, 116, 112, 110, 108, 106, + 123, 126, 129, 131, 133, 135, 137, 138, 139, 141, 143, 147, 150, 156, 161, 164, 167, 170, 173, 177, 181, 184, 187, 188, 190, 189, 187, 185, 183, 179, 176, 174, 174, 174, 174, 174, 176, 177, 179, 180, 182, 183, 182, 181, 181, 180, 176, 171, 166, 160, 152, 147, 142, 138, 133, 126, 121, 115, 110, 106, 105, + 124, 127, 130, 132, 135, 137, 138, 140, 142, 144, 147, 149, 154, 157, 161, 165, 168, 171, 175, 178, 181, 184, 186, 187, 187, 184, 184, 181, 179, 175, 171, 169, 168, 168, 168, 169, 170, 172, 174, 177, 178, 179, 180, 181, 181, 180, 179, 174, 167, 161, 155, 148, 144, 139, 134, 128, 121, 115, 110, 106, 105, + 123, 128, 131, 133, 136, 138, 140, 142, 144, 146, 149, 151, 154, 157, 160, 164, 168, 172, 175, 178, 181, 183, 184, 184, 185, 183, 180, 177, 174, 170, 167, 165, 164, 164, 164, 165, 166, 168, 171, 175, 176, 178, 180, 181, 180, 180, 179, 177, 170, 163, 157, 150, 144, 139, 134, 128, 121, 115, 110, 108, 107, + 123, 127, 131, 134, 136, 138, 140, 142, 144, 147, 149, 151, 154, 157, 160, 164, 168, 171, 174, 178, 180, 181, 181, 182, 183, 181, 178, 173, 169, 166, 163, 161, 161, 160, 160, 161, 163, 165, 168, 173, 176, 178, 179, 180, 181, 180, 180, 175, 173, 166, 159, 152, 145, 139, 134, 127, 121, 115, 110, 109, 108, + 120, 124, 128, 131, 134, 137, 139, 142, 144, 146, 149, 151, 153, 156, 160, 163, 167, 171, 174, 178, 180, 180, 180, 180, 180, 180, 175, 171, 167, 162, 160, 158, 157, 157, 157, 158, 159, 162, 166, 170, 175, 177, 178, 180, 181, 181, 180, 178, 175, 169, 160, 154, 148, 140, 134, 128, 121, 115, 110, 110, 109, + 118, 121, 125, 129, 132, 134, 137, 140, 142, 145, 147, 149, 151, 155, 159, 163, 166, 169, 173, 177, 179, 180, 180, 180, 180, 179, 174, 169, 166, 161, 158, 156, 154, 153, 153, 154, 156, 159, 163, 169, 173, 175, 178, 180, 181, 180, 180, 179, 175, 170, 160, 154, 149, 142, 135, 128, 122, 116, 111, 110, 110, + 117, 120, 121, 125, 129, 132, 135, 138, 140, 143, 145, 147, 149, 153, 157, 160, 163, 166, 171, 174, 177, 179, 180, 180, 180, 179, 172, 168, 164, 160, 157, 154, 151, 149, 150, 150, 154, 158, 164, 169, 174, 178, 180, 180, 180, 180, 178, 177, 175, 170, 161, 153, 148, 142, 135, 129, 123, 116, 113, 112, 110, + 115, 118, 120, 122, 126, 130, 133, 136, 138, 141, 143, 145, 148, 151, 154, 157, 160, 163, 168, 171, 174, 177, 179, 179, 179, 176, 171, 167, 164, 160, 156, 153, 149, 148, 149, 151, 155, 158, 163, 170, 173, 177, 179, 180, 180, 180, 178, 175, 173, 171, 162, 154, 147, 141, 136, 130, 124, 117, 115, 112, 110, + 114, 116, 118, 120, 122, 127, 131, 133, 136, 138, 141, 143, 146, 148, 151, 154, 157, 160, 164, 168, 171, 174, 178, 178, 179, 177, 173, 169, 165, 161, 157, 154, 151, 149, 150, 152, 155, 159, 166, 171, 175, 177, 179, 180, 180, 179, 176, 174, 171, 168, 159, 151, 146, 141, 135, 129, 124, 119, 116, 113, 110, + 115, 114, 116, 118, 120, 122, 127, 129, 132, 136, 139, 141, 143, 146, 148, 151, 153, 156, 160, 164, 167, 172, 174, 176, 177, 176, 173, 170, 166, 162, 159, 157, 154, 153, 154, 155, 158, 161, 169, 172, 174, 176, 178, 178, 178, 178, 175, 172, 169, 162, 156, 149, 144, 140, 134, 128, 123, 118, 115, 112, 110, + 113, 113, 114, 116, 118, 120, 122, 125, 129, 133, 136, 138, 141, 143, 146, 149, 150, 153, 156, 160, 165, 170, 173, 176, 176, 176, 173, 172, 169, 165, 163, 160, 158, 157, 158, 159, 161, 166, 170, 170, 173, 175, 176, 178, 176, 173, 171, 168, 164, 158, 153, 146, 140, 137, 132, 127, 121, 117, 113, 111, 110, + 111, 112, 113, 114, 116, 118, 120, 122, 126, 130, 133, 136, 139, 142, 145, 147, 148, 151, 155, 158, 163, 168, 173, 176, 177, 177, 176, 174, 171, 169, 166, 164, 161, 161, 162, 164, 165, 167, 170, 170, 171, 173, 173, 173, 170, 168, 165, 163, 160, 155, 149, 143, 138, 134, 130, 125, 119, 116, 112, 110, 109, + 110, 112, 113, 113, 114, 116, 118, 120, 123, 127, 131, 134, 137, 141, 143, 145, 148, 150, 154, 157, 161, 166, 171, 176, 178, 178, 178, 176, 174, 172, 170, 167, 167, 167, 166, 168, 170, 169, 168, 167, 168, 168, 168, 168, 167, 165, 163, 160, 156, 152, 146, 140, 136, 131, 128, 122, 118, 114, 110, 110, 109, + 109, 110, 111, 112, 114, 116, 118, 119, 120, 124, 128, 131, 136, 140, 142, 145, 147, 150, 153, 157, 160, 165, 170, 174, 178, 179, 179, 178, 178, 176, 174, 171, 170, 170, 170, 168, 167, 166, 164, 163, 161, 162, 163, 163, 163, 161, 160, 157, 153, 148, 142, 136, 130, 127, 124, 120, 117, 113, 110, 110, 109, + 108, 109, 111, 112, 114, 116, 117, 118, 120, 121, 125, 128, 132, 138, 142, 144, 147, 149, 153, 156, 160, 164, 170, 174, 178, 180, 180, 179, 179, 178, 176, 172, 170, 170, 170, 168, 166, 164, 162, 160, 157, 156, 157, 158, 158, 156, 153, 151, 149, 144, 139, 130, 127, 124, 121, 118, 115, 112, 110, 110, 109, + 108, 109, 111, 113, 114, 116, 117, 118, 119, 120, 122, 126, 130, 135, 139, 143, 147, 149, 152, 156, 160, 164, 169, 173, 177, 180, 180, 180, 180, 179, 178, 174, 170, 170, 168, 167, 165, 163, 161, 157, 154, 153, 152, 152, 152, 149, 148, 147, 144, 140, 134, 128, 125, 122, 119, 117, 114, 110, 110, 109, 109, + 107, 108, 111, 112, 114, 115, 116, 117, 119, 120, 121, 124, 128, 133, 137, 141, 145, 149, 152, 156, 160, 164, 168, 172, 176, 179, 180, 180, 180, 179, 178, 174, 170, 168, 166, 165, 163, 161, 158, 154, 150, 149, 148, 146, 145, 143, 143, 143, 140, 136, 130, 126, 123, 120, 118, 115, 112, 110, 110, 109, 109, + 107, 108, 110, 112, 113, 113, 115, 116, 118, 120, 122, 125, 128, 132, 136, 140, 145, 148, 150, 155, 160, 164, 167, 170, 174, 177, 179, 179, 178, 176, 176, 173, 169, 166, 164, 163, 161, 159, 155, 152, 148, 145, 143, 141, 140, 139, 139, 138, 136, 132, 128, 124, 121, 118, 116, 114, 111, 110, 110, 109, 108, + 107, 108, 109, 111, 113, 114, 116, 117, 119, 120, 122, 125, 128, 132, 137, 141, 144, 146, 149, 152, 157, 162, 166, 168, 171, 173, 175, 175, 173, 172, 172, 171, 168, 165, 162, 160, 158, 156, 153, 149, 145, 142, 139, 138, 137, 136, 135, 133, 131, 129, 126, 122, 119, 117, 114, 112, 110, 110, 109, 108, 107, + 108, 109, 110, 112, 114, 115, 116, 117, 119, 120, 122, 126, 129, 133, 137, 141, 143, 146, 148, 151, 155, 160, 164, 167, 168, 169, 170, 170, 169, 168, 167, 168, 166, 163, 160, 158, 155, 153, 150, 147, 143, 140, 137, 136, 134, 133, 132, 130, 129, 127, 125, 121, 118, 115, 112, 110, 110, 110, 108, 107, 107, + 109, 110, 111, 113, 115, 116, 117, 118, 120, 121, 123, 126, 129, 133, 138, 141, 143, 146, 148, 150, 155, 159, 163, 165, 166, 167, 168, 168, 166, 165, 164, 161, 160, 159, 158, 155, 152, 149, 147, 144, 141, 138, 135, 134, 132, 130, 129, 128, 126, 124, 122, 120, 117, 113, 111, 110, 110, 110, 108, 107, 107, + 110, 111, 112, 113, 116, 117, 118, 119, 120, 122, 125, 127, 130, 133, 138, 141, 143, 146, 148, 150, 154, 159, 162, 163, 164, 166, 166, 166, 165, 163, 161, 159, 157, 156, 155, 153, 150, 146, 143, 140, 138, 136, 133, 132, 130, 129, 128, 125, 124, 122, 120, 119, 117, 114, 111, 110, 110, 109, 108, 107, 107, + 111, 112, 113, 114, 116, 117, 118, 119, 120, 123, 125, 128, 130, 134, 139, 141, 144, 146, 148, 151, 154, 158, 161, 164, 166, 167, 168, 166, 165, 163, 161, 158, 156, 154, 152, 150, 146, 142, 139, 137, 135, 133, 131, 130, 129, 128, 127, 125, 123, 121, 120, 118, 116, 113, 111, 110, 110, 109, 108, 107, 106, + 111, 112, 113, 115, 117, 118, 118, 120, 121, 124, 126, 128, 131, 135, 139, 142, 144, 146, 148, 151, 155, 160, 164, 165, 168, 169, 169, 168, 166, 163, 160, 158, 156, 153, 151, 148, 145, 142, 139, 137, 135, 132, 130, 129, 127, 126, 125, 124, 123, 120, 120, 117, 116, 114, 112, 110, 110, 108, 107, 106, 106, + 112, 113, 114, 116, 117, 118, 119, 120, 122, 124, 127, 129, 132, 135, 139, 142, 144, 146, 149, 152, 157, 162, 167, 169, 170, 170, 170, 168, 165, 163, 161, 159, 157, 155, 151, 148, 145, 141, 139, 136, 134, 132, 130, 128, 127, 126, 124, 123, 122, 120, 119, 117, 116, 114, 112, 111, 109, 107, 106, 106, 105, + 113, 114, 115, 116, 117, 119, 119, 120, 122, 125, 127, 129, 132, 135, 139, 142, 144, 147, 149, 154, 159, 164, 169, 170, 170, 170, 170, 170, 168, 165, 163, 161, 158, 155, 151, 148, 145, 142, 139, 137, 135, 132, 131, 128, 126, 125, 124, 122, 121, 120, 119, 117, 115, 113, 111, 110, 109, 106, 105, 105, 104, + 113, 114, 115, 117, 118, 119, 120, 121, 123, 125, 127, 130, 132, 135, 139, 142, 145, 148, 150, 156, 161, 166, 170, 170, 170, 170, 170, 170, 169, 166, 163, 161, 159, 155, 151, 148, 146, 143, 140, 138, 135, 134, 132, 130, 127, 125, 123, 121, 120, 120, 119, 116, 114, 112, 110, 110, 108, 106, 105, 104, 104, + 114, 115, 116, 117, 118, 119, 120, 121, 123, 126, 128, 130, 133, 136, 139, 142, 145, 148, 152, 157, 161, 166, 168, 170, 170, 170, 170, 168, 166, 164, 163, 160, 159, 155, 151, 148, 146, 143, 141, 138, 136, 134, 132, 130, 128, 125, 123, 121, 120, 120, 118, 116, 113, 111, 110, 110, 109, 106, 105, 104, 104, + 115, 116, 117, 118, 119, 120, 121, 121, 123, 126, 128, 131, 134, 136, 139, 142, 145, 149, 152, 157, 161, 163, 164, 166, 168, 167, 166, 164, 163, 161, 160, 158, 156, 152, 149, 147, 144, 143, 141, 139, 136, 134, 132, 130, 128, 125, 122, 120, 120, 119, 117, 115, 113, 110, 110, 109, 107, 106, 105, 104, 104, + 115, 116, 117, 118, 119, 120, 121, 122, 123, 125, 128, 131, 134, 137, 139, 142, 145, 149, 152, 156, 159, 159, 160, 162, 162, 161, 161, 160, 159, 158, 157, 155, 153, 150, 148, 146, 145, 143, 142, 140, 137, 134, 131, 129, 126, 124, 122, 120, 119, 117, 115, 113, 111, 110, 109, 109, 107, 106, 105, 104, 104, + 114, 115, 116, 116, 118, 119, 120, 121, 122, 126, 129, 132, 135, 137, 140, 143, 146, 149, 152, 155, 156, 157, 158, 159, 159, 159, 158, 158, 157, 155, 153, 151, 150, 149, 147, 146, 145, 144, 142, 141, 138, 135, 132, 128, 125, 122, 120, 118, 117, 115, 113, 112, 110, 109, 108, 108, 106, 105, 105, 104, 104, + 113, 114, 115, 116, 117, 118, 119, 120, 123, 126, 129, 132, 135, 138, 140, 143, 146, 148, 151, 153, 154, 156, 157, 157, 157, 157, 156, 155, 154, 152, 150, 149, 148, 147, 146, 145, 144, 142, 141, 140, 139, 136, 132, 129, 125, 121, 118, 116, 115, 113, 111, 110, 109, 108, 108, 107, 106, 105, 104, 104, 104, + 112, 113, 114, 115, 116, 117, 119, 120, 122, 126, 130, 133, 136, 138, 141, 143, 146, 148, 150, 152, 154, 155, 155, 155, 155, 155, 154, 152, 152, 150, 148, 147, 146, 145, 145, 143, 142, 141, 140, 140, 140, 137, 133, 129, 125, 120, 117, 115, 111, 110, 110, 109, 108, 107, 107, 106, 105, 105, 104, 104, 103, + 111, 112, 114, 115, 116, 117, 118, 120, 122, 125, 131, 134, 137, 139, 142, 144, 146, 148, 150, 152, 153, 153, 153, 153, 153, 153, 153, 151, 149, 147, 146, 144, 144, 143, 143, 142, 141, 140, 140, 140, 140, 138, 134, 130, 123, 120, 118, 111, 110, 110, 110, 108, 107, 106, 108, 105, 105, 104, 104, 103, 103, + 111, 112, 113, 115, 115, 116, 117, 119, 121, 126, 131, 135, 138, 140, 142, 144, 146, 148, 150, 151, 151, 151, 151, 151, 151, 151, 151, 150, 148, 146, 144, 142, 141, 141, 142, 141, 140, 140, 140, 140, 140, 140, 136, 132, 126, 120, 115, 110, 110, 110, 109, 107, 106, 105, 107, 105, 104, 104, 104, 103, 103, + 112, 113, 113, 114, 115, 116, 117, 119, 122, 127, 132, 135, 139, 141, 143, 145, 147, 149, 150, 150, 150, 150, 150, 150, 150, 150, 150, 149, 147, 144, 142, 141, 140, 140, 140, 140, 140, 140, 140, 140, 140, 140, 137, 133, 128, 120, 117, 110, 110, 110, 108, 106, 105, 105, 106, 105, 104, 104, 103, 103, 103, + 112, 113, 114, 114, 116, 117, 118, 120, 122, 128, 132, 136, 139, 141, 144, 146, 147, 149, 150, 150, 150, 150, 150, 150, 150, 150, 150, 149, 146, 143, 141, 140, 140, 139, 139, 139, 140, 140, 140, 140, 140, 140, 137, 133, 129, 121, 118, 110, 110, 109, 107, 106, 105, 105, 105, 104, 104, 103, 103, 103, 102, + 112, 114, 114, 115, 116, 117, 119, 120, 122, 128, 133, 136, 140, 142, 144, 146, 148, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 148, 145, 142, 140, 138, 138, 138, 137, 138, 140, 140, 140, 140, 140, 140, 137, 134, 130, 122, 118, 110, 110, 108, 106, 105, 103, 104, 104, 104, 104, 103, 103, 102, 102, + 113, 114, 115, 116, 116, 117, 118, 120, 123, 129, 133, 137, 140, 142, 144, 146, 149, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 147, 143, 141, 139, 137, 136, 136, 135, 136, 138, 140, 140, 140, 140, 139, 136, 134, 130, 123, 119, 113, 109, 108, 106, 104, 103, 104, 104, 104, 103, 103, 102, 102, 101, + 114, 115, 115, 116, 117, 118, 118, 120, 123, 129, 133, 137, 140, 143, 145, 147, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 148, 145, 142, 139, 138, 136, 135, 134, 134, 134, 136, 138, 137, 138, 139, 137, 134, 132, 125, 122, 117, 114, 109, 107, 105, 103, 102, 104, 104, 103, 103, 102, 102, 101, 101, + 114, 115, 116, 117, 117, 119, 118, 120, 123, 128, 132, 136, 139, 142, 145, 148, 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 147, 144, 141, 139, 136, 135, 134, 133, 132, 132, 134, 134, 134, 134, 135, 133, 131, 128, 124, 120, 116, 113, 110, 107, 104, 102, 102, 103, 103, 103, 102, 102, 102, 101, 100, + 115, 116, 116, 117, 118, 119, 119, 120, 124, 128, 132, 136, 139, 142, 145, 148, 150, 150, 150, 150, 150, 150, 150, 150, 150, 149, 146, 143, 140, 138, 135, 134, 133, 131, 131, 131, 131, 131, 131, 131, 130, 127, 124, 122, 119, 117, 115, 112, 109, 106, 104, 101, 102, 103, 103, 102, 102, 102, 101, 100, 100, + 115, 116, 117, 118, 118, 119, 120, 123, 125, 128, 131, 135, 138, 141, 145, 148, 150, 150, 150, 150, 150, 150, 150, 150, 150, 147, 145, 142, 139, 137, 134, 132, 131, 130, 129, 128, 128, 128, 128, 128, 126, 123, 121, 119, 116, 114, 112, 110, 108, 105, 103, 101, 103, 103, 103, 102, 102, 101, 100, 100, 100, + 116, 117, 118, 118, 119, 120, 122, 123, 125, 128, 131, 134, 137, 141, 145, 148, 149, 150, 150, 150, 150, 150, 150, 150, 148, 145, 143, 141, 138, 135, 133, 130, 129, 128, 127, 126, 125, 125, 125, 124, 123, 120, 118, 116, 114, 111, 109, 107, 106, 104, 102, 100, 101, 101, 102, 102, 101, 100, 100, 100, 100, + 116, 117, 118, 119, 120, 121, 123, 124, 126, 128, 130, 133, 137, 140, 144, 145, 147, 148, 149, 150, 149, 149, 147, 146, 144, 141, 139, 136, 133, 131, 129, 128, 127, 126, 125, 124, 123, 123, 122, 121, 120, 118, 116, 114, 112, 108, 107, 105, 103, 102, 100, 100, 100, 100, 101, 101, 100, 100, 100, 100, 100, + 117, 118, 119, 119, 120, 121, 123, 124, 126, 128, 129, 131, 135, 139, 142, 143, 145, 146, 147, 147, 147, 146, 144, 142, 140, 138, 135, 133, 130, 128, 127, 126, 125, 124, 123, 122, 121, 120, 119, 118, 117, 115, 114, 112, 110, 106, 105, 102, 101, 100, 100, 100, 100, 100, 100, 100, 100, 99, 99, 99, 99, + 117, 118, 119, 120, 120, 121, 123, 124, 125, 126, 128, 129, 132, 137, 140, 142, 143, 143, 144, 144, 144, 143, 141, 139, 137, 135, 133, 130, 128, 127, 126, 125, 123, 122, 121, 120, 119, 117, 116, 115, 114, 112, 111, 108, 107, 105, 100, 100, 100, 100, 100, 100, 100, 99, 99, 99, 99, 99, 99, 99, 98, + 116, 117, 118, 120, 120, 121, 122, 123, 124, 125, 126, 128, 130, 134, 139, 140, 141, 141, 141, 141, 141, 140, 138, 136, 134, 133, 131, 129, 127, 125, 124, 123, 122, 120, 119, 118, 117, 116, 114, 112, 111, 108, 109, 106, 106, 100, 100, 100, 100, 100, 99, 99, 99, 99, 99, 99, 99, 98, 98, 98, 97, + 114, 115, 116, 117, 119, 119, 120, 121, 122, 123, 125, 127, 129, 133, 136, 134, 134, 136, 138, 138, 137, 137, 135, 133, 132, 130, 129, 127, 125, 124, 122, 121, 120, 119, 117, 116, 115, 114, 112, 110, 109, 108, 107, 105, 105, 100, 100, 100, 100, 99, 99, 99, 98, 98, 98, 98, 98, 97, 97, 97, 97, + 112, 113, 114, 115, 116, 116, 117, 119, 120, 122, 124, 126, 127, 129, 129, 128, 127, 129, 132, 133, 133, 133, 133, 131, 129, 127, 126, 125, 124, 122, 121, 119, 118, 117, 116, 114, 113, 112, 110, 109, 108, 106, 106, 105, 100, 100, 100, 98, 98, 98, 98, 98, 98, 97, 97, 97, 97, 97, 97, 97, 96, + 109, 111, 112, 112, 113, 113, 113, 114, 116, 119, 121, 123, 124, 125, 124, 123, 123, 123, 125, 127, 129, 129, 128, 128, 127, 125, 124, 123, 122, 121, 119, 118, 117, 116, 114, 113, 112, 110, 109, 108, 107, 106, 105, 100, 100, 100, 97, 97, 97, 97, 97, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, + 106, 107, 108, 108, 109, 110, 110, 112, 113, 114, 117, 119, 120, 121, 119, 117, 117, 117, 118, 120, 123, 124, 125, 125, 125, 123, 121, 120, 120, 119, 118, 117, 116, 115, 114, 113, 111, 109, 109, 107, 106, 105, 100, 100, 100, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, + 104, 105, 105, 106, 106, 107, 108, 108, 109, 109, 111, 115, 116, 114, 113, 112, 111, 110, 111, 113, 116, 119, 122, 122, 122, 121, 120, 119, 118, 118, 117, 116, 115, 114, 113, 112, 111, 108, 108, 106, 105, 100, 100, 100, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, + 102, 103, 103, 104, 104, 105, 106, 106, 107, 108, 109, 111, 112, 110, 109, 108, 108, 108, 108, 109, 110, 112, 116, 117, 117, 118, 118, 118, 117, 116, 116, 115, 114, 113, 112, 111, 110, 107, 107, 105, 100, 100, 100, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, + 101, 102, 103, 103, 104, 105, 105, 106, 106, 107, 108, 109, 109, 107, 106, 106, 105, 105, 105, 106, 107, 108, 109, 110, 111, 113, 114, 115, 115, 115, 114, 113, 112, 111, 110, 108, 108, 106, 105, 100, 100, 100, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, + 100, 101, 102, 102, 103, 103, 104, 104, 105, 106, 106, 107, 106, 106, 106, 105, 105, 104, 103, 103, 104, 105, 107, 108, 110, 111, 111, 112, 112, 113, 113, 112, 111, 110, 108, 107, 106, 105, 100, 100, 100, 98, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, + 100, 101, 101, 102, 102, 103, 103, 104, 104, 105, 105, 105, 105, 106, 105, 105, 104, 103, 102, 101, 102, 103, 104, 106, 107, 110, 111, 111, 111, 112, 112, 112, 110, 107, 107, 106, 105, 102, 100, 100, 99, 98, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 95, + 99, 100, 101, 102, 102, 103, 103, 103, 104, 104, 104, 104, 103, 104, 104, 104, 104, 102, 101, 101, 102, 103, 104, 105, 107, 110, 111, 111, 111, 111, 111, 111, 108, 106, 105, 105, 102, 101, 100, 99, 99, 98, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 95, 95, + 99, 100, 100, 101, 101, 102, 102, 102, 103, 103, 103, 103, 102, 103, 103, 104, 103, 102, 101, 101, 101, 102, 103, 104, 106, 109, 110, 111, 111, 111, 110, 110, 107, 105, 103, 104, 100, 100, 99, 99, 98, 98, 97, 97, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 95, 95, 95, 95, 95, 95, 95, + 99, 100, 100, 100, 101, 101, 101, 102, 102, 103, 102, 102, 101, 102, 102, 103, 103, 101, 101, 100, 101, 101, 102, 103, 105, 109, 110, 110, 111, 110, 110, 109, 106, 105, 100, 102, 100, 99, 99, 99, 98, 98, 97, 97, 96, 96, 96, 96, 96, 96, 95, 95, 95, 95, 95, 95, 95, 95, 95, 95, 94, + 99, 99, 99, 99, 100, 100, 101, 101, 102, 102, 101, 101, 101, 101, 101, 102, 102, 101, 100, 100, 101, 101, 101, 103, 104, 107, 109, 109, 110, 110, 109, 108, 105, 102, 100, 100, 99, 99, 99, 98, 98, 98, 97, 96, 96, 96, 96, 96, 95, 95, 95, 95, 95, 95, 95, 94, 94, 94, 94, 94, 94, + 98, 99, 99, 99, 99, 100, 100, 101, 101, 102, 101, 100, 100, 100, 101, 101, 101, 100, 100, 100, 100, 101, 101, 101, 103, 106, 107, 109, 109, 109, 109, 107, 104, 101, 100, 99, 99, 99, 98, 98, 98, 97, 96, 96, 96, 96, 95, 95, 95, 95, 95, 95, 95, 94, 94, 94, 94, 94, 94, 94, 94, + 98, 98, 98, 99, 99, 99, 100, 100, 101, 101, 100, 100, 99, 99, 100, 100, 100, 100, 100, 100, 100, 101, 101, 101, 102, 105, 106, 109, 108, 109, 107, 105, 102, 100, 100, 99, 99, 98, 98, 98, 97, 96, 96, 96, 96, 95, 95, 95, 95, 95, 95, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, + 97, 98, 98, 98, 99, 99, 99, 100, 100, 100, 100, 100, 99, 99, 99, 100, 100, 100, 100, 100, 100, 100, 101, 101, 101, 103, 104, 105, 106, 105, 104, 101, 100, 100, 99, 99, 98, 98, 97, 97, 97, 96, 96, 96, 95, 95, 95, 95, 95, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, + 97, 97, 97, 98, 98, 99, 99, 99, 100, 100, 100, 99, 99, 99, 99, 99, 100, 100, 100, 100, 100, 100, 101, 101, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 99, 99, 98, 97, 97, 97, 96, 96, 96, 95, 95, 95, 95, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, 94, +})}