Files
tensor/plot/plot.go
2026-09-28 16:52:31 +02:00

273 lines
9.3 KiB
Go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT
// Package plot draws the deterministic SVG line charts a scientific
// paper needs: linear axes, five ticks each, one legend line per
// series, and nothing else. The output is deterministic by contract:
// the same chart always renders byte for byte the same file, so a
// figure in a paper can be regenerated and compared exactly like any
// other computed number. The package is small by intent; it draws the
// figures, it does not stage a cinema.
package plot
import (
"fmt"
"math"
"os"
"path/filepath"
"strings"
"sourcedock.dev/petrbalvin/tensor/internal/base"
"sourcedock.dev/petrbalvin/tensor/internal/core"
)
// Point is one data point in axis units.
type Point struct {
X, Y float64
}
// Series is one named polyline.
type Series struct {
Name string
Points []Point
}
// Chart is a linear-axis line chart.
type Chart struct {
Title string
XLabel string
YLabel string
Width int
Height int
Series []Series
// XRange and YRange are optional; a zero, inverted or non-finite
// span falls back to the data's own bounds.
XRange [2]float64
YRange [2]float64
}
// Line returns a series joining the points (xs[i], ys[i]) of two
// arrays. Both arrays must be rank 1, of equal, non-zero length, and
// hold only finite numbers; the values are read through the promotion
// ladder, so any numeric dtype is accepted.
func Line(name string, xs, ys *core.Array) (Series, error) {
const op = "Plot"
if xs == nil || ys == nil {
return Series{}, base.Errf("%s: a nil array cannot make a series", op)
}
if xs.NDim() != 1 {
return Series{}, base.Errf("%s: the x values must be rank 1, got shape %s", op, base.ShapeText(xs.Shape()))
}
if ys.NDim() != 1 {
return Series{}, base.Errf("%s: the y values must be rank 1, got shape %s", op, base.ShapeText(ys.Shape()))
}
if xs.Len() != ys.Len() {
return Series{}, base.Errf("%s: length mismatch, %d points of x against %d points of y",
op, xs.Len(), ys.Len())
}
if xs.Len() == 0 {
return Series{}, base.Errf("%s: an empty array cannot make a series", op)
}
pts := make([]Point, xs.Len())
for i := range pts {
x, y := xs.FloatAt(i), ys.FloatAt(i)
if math.IsNaN(x) || math.IsInf(x, 0) || math.IsNaN(y) || math.IsInf(y, 0) {
return Series{}, base.Errf("%s: non-finite point at index %d: (%g, %g)", op, i, x, y)
}
pts[i] = Point{X: x, Y: y}
}
return Series{Name: name, Points: pts}, nil
}
// WriteSVG renders the chart into path. The output is deterministic:
// the same chart always renders byte for byte the same file. Every
// point of every series must be finite, the contract the Line
// constructor enforces on the caller's behalf and this entry point
// enforces for a series built by hand.
func (c Chart) WriteSVG(path string) error {
w := c.Width
if w <= 0 {
w = 720
}
h := c.Height
if h <= 0 {
h = 460
}
const (
left = 64.0
right = 16.0
top = 40.0
bottom = 52.0
)
all := make([]Point, 0, 256)
for si, s := range c.Series {
for pi, p := range s.Points {
if math.IsNaN(p.X) || math.IsInf(p.X, 0) || math.IsNaN(p.Y) || math.IsInf(p.Y, 0) {
return base.Errf("Plot: series %d (%s) holds the non-finite point %d: (%g, %g)", si, s.Name, pi, p.X, p.Y)
}
}
all = append(all, s.Points...)
}
if len(all) < 2 {
return base.Errf("Plot: the chart needs at least two points, has %d", len(all))
}
xr := c.XRange
if !(xr[0] < xr[1]) || math.IsInf(xr[0], 0) || math.IsInf(xr[1], 0) {
xr = bounds(all, true)
}
yr := c.YRange
if !(yr[0] < yr[1]) || math.IsInf(yr[0], 0) || math.IsInf(yr[1], 0) {
yr = bounds(all, false)
}
// project maps a value of the range [r0, r1] onto the plot span
// starting at plotLo. The normal form is the arithmetic the chart
// has always run; a span, or an offset from the range's start, that
// overflows float64 falls back to the halved form, whose every term
// a finite input keeps finite, so a range the data spans but the
// arithmetic cannot still draws instead of carrying a coordinate no
// renderer displays.
project := func(v, r0, r1, plotLo, plotSpan float64) float64 {
u := (v - r0) / (r1 - r0)
if math.IsNaN(u) || math.IsInf(u, 0) {
u = (v/2 - r0/2) / (r1/2 - r0/2)
}
return plotLo + u*plotSpan
}
px := func(x float64) float64 {
return project(x, xr[0], xr[1], left, float64(w)-left-right)
}
py := func(y float64) float64 {
return float64(h) - bottom - project(y, yr[0], yr[1], 0, float64(h)-top-bottom)
}
// tickValue places the k-th of the five ticks. A span that overflows
// float64 makes the affine form NaN or Inf, so the tick falls back to
// the convex combination, which stays between the range's own finite
// ends.
tickValue := func(r0, r1 float64, k int) float64 {
t := r0 + (r1-r0)*float64(k)/4
if math.IsNaN(t) || math.IsInf(t, 0) {
f := float64(k) / 4
t = r0*(1-f) + r1*f
}
return t
}
var b strings.Builder
b.WriteString(xmlHeader)
fmt.Fprintf(&b, "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" viewBox=\"0 0 %d %d\">\n", w, h, w, h)
fmt.Fprintf(&b, "<rect width=\"%d\" height=\"%d\" fill=\"white\"/>\n", w, h)
fmt.Fprintf(&b, "<text x=\"%g\" y=\"24\" font-family=\"sans-serif\" font-size=\"15\" fill=\"#111\">%s</text>\n",
left, esc(c.Title))
// Axes with five ticks each.
for k := range 5 {
t := tickValue(xr[0], xr[1], k)
x := px(t)
fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#ccc\" stroke-width=\"1\"/>\n",
x, top, x, float64(h)-bottom)
fmt.Fprintf(&b, "<text x=\"%g\" y=\"%g\" font-family=\"sans-serif\" font-size=\"11\" fill=\"#333\" text-anchor=\"middle\">%s</text>\n",
x, float64(h)-bottom+16, tick(t))
}
for k := range 5 {
t := tickValue(yr[0], yr[1], k)
y := py(t)
fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#ccc\" stroke-width=\"1\"/>\n",
left, y, float64(w)-right, y)
fmt.Fprintf(&b, "<text x=\"%g\" y=\"%g\" font-family=\"sans-serif\" font-size=\"11\" fill=\"#333\" text-anchor=\"end\">%s</text>\n",
left-6, y+4, tick(t))
}
fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#111\" stroke-width=\"1\"/>\n",
left, float64(h)-bottom, float64(w)-right, float64(h)-bottom)
fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"#111\" stroke-width=\"1\"/>\n",
left, top, left, float64(h)-bottom)
fmt.Fprintf(&b, "<text x=\"%g\" y=\"%g\" font-family=\"sans-serif\" font-size=\"12\" fill=\"#111\" text-anchor=\"middle\">%s</text>\n",
(left+float64(w)-right)/2, float64(h)-12, esc(c.XLabel))
fmt.Fprintf(&b, "<text x=\"16\" y=\"%g\" font-family=\"sans-serif\" font-size=\"12\" fill=\"#111\">%s</text>\n",
top-12, esc(c.YLabel))
for i, s := range c.Series {
colour := colour(i)
fmt.Fprintf(&b, "<polyline fill=\"none\" stroke=\"%s\" stroke-width=\"1.8\" points=\"", colour)
for j, p := range s.Points {
if j > 0 {
b.WriteByte(' ')
}
fmt.Fprintf(&b, "%.2f,%.2f", px(p.X), py(p.Y))
}
b.WriteString("\"/>\n")
ly := top + 16 + float64(i)*16
fmt.Fprintf(&b, "<line x1=\"%g\" y1=\"%g\" x2=\"%g\" y2=\"%g\" stroke=\"%s\" stroke-width=\"1.8\"/>\n",
float64(w)-230, ly, float64(w)-214, ly, colour)
fmt.Fprintf(&b, "<text x=\"%g\" y=\"%g\" font-family=\"sans-serif\" font-size=\"11\" fill=\"#111\">%s</text>\n",
float64(w)-208, ly+4, esc(s.Name))
}
b.WriteString("</svg>\n")
if err := os.MkdirAll(filepath.Dir(path), 0o755); err != nil {
return base.Errf("Plot: %w", err)
}
if err := os.WriteFile(path, []byte(b.String()), 0o644); err != nil {
return base.Errf("Plot: %w", err)
}
return nil
}
func bounds(pts []Point, xAxis bool) [2]float64 {
lo, hi := math.Inf(1), math.Inf(-1)
for _, p := range pts {
v := p.Y
if xAxis {
v = p.X
}
lo = math.Min(lo, v)
hi = math.Max(hi, v)
}
if hi == lo {
hi = lo + 1
}
padLo, padHi := lo-0.05*(hi-lo), hi+0.05*(hi-lo)
if math.IsInf(padLo, 0) || math.IsInf(padHi, 0) {
// The padding, or the span it scales, overflows the range the
// finite data itself fits; the unpadded bounds keep every
// projection finite.
return [2]float64{lo, hi}
}
return [2]float64{padLo, padHi}
}
func tick(v float64) string {
if v == math.Trunc(v) && math.Abs(v) < 1e15 {
return fmt.Sprintf("%d", int64(v))
}
return fmt.Sprintf("%g", v)
}
func esc(s string) string {
r := strings.NewReplacer("&", "&amp;", "<", "&lt;", ">", "&gt;", `"`, "&#34;")
return r.Replace(s)
}
// seriesColours is the chart's fixed colour cycle: seven even samples
// of the Viridis perceptual-uniform map (Nathaniel J. Smith, Stéfan
// van der Walt and Eric Firing, released under CC0), read from the
// map's 256-entry table at t = 0, 7/60, ..., 0.7 by linear
// interpolation, each channel rounded to the nearest byte. The map's
// light tail is left out on purpose: the
// chart paints on white, and the pale yellows the full range ends in
// drop far below a legible contrast at stroke width, while the
// sampled range runs dark violet through blue and teal to green with
// every stroke legible. The cycle is a constant, so the same chart
// renders the same colours byte for byte, like everything else it
// draws.
var seriesColours = [7]string{
"#440154",
"#482a79",
"#3d4d8a",
"#2f6c8e",
"#23888e",
"#20a486",
"#43bf71",
}
func colour(i int) string {
return seriesColours[i%len(seriesColours)]
}
const xmlHeader = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"