nnls.go

package main

import (
	"fmt"
	"log"
	"math"
	"math/rand"
	"os"
	"testing"
	"time"

	"github.com/gonum/plot"
	"github.com/gonum/plot/plotter"
	"github.com/gonum/plot/plotutil"
	"github.com/gonum/plot/vg"
	reg "github.com/sajari/regression"
	"github.com/soniakeys/graph"
	"github.com/soniakeys/graph/dot"
	"github.com/soniakeys/nnls"
)

func main() {
	var l reg.Regression
	l.SetObservedName("ns/op")
	//	l.SetVarName(0, "n")
	//	l.SetVarName(1, "log n")
	//	l.SetVarName(2, "log^2 n")
	l.SetVarName(0, "n log n")
	//	l.SetVarName(4, "n^2")
	r := rand.New(rand.NewSource(time.Now().Unix()))
	ns := []int{16, 32, 64, 128, 256, 512, 1024, 2048, 4096,
		8192, 16384, 32768, 65536, 131072, 262144, 524288, 1048576,
		2097152, 4194304} /*, 8388608, 16777216}*/
	rep := 5
	for _, n := range ns {
		g, pos, wt, err := graph.LabeledEuclidean(n, n*10, 4, 100, r)
		if err != nil {
			log.Fatal(err)
		}

		{
			f, err := os.Create(fmt.Sprintf("bench%d.dot", n))
			if err != nil {
				log.Fatal(err)
			}
			dot.Write(g, f, dot.NodePos(func(n graph.NI) string {
				return fmt.Sprintf("%.2f,%.2f", 10*pos[n].X, 10*pos[n].Y)
			}))
			f.Close()
		}
		d := graph.NewDijkstra(g.LabeledAdjacencyList, func(l graph.LI) float64 {
			return wt[l]
		})
		for j := 0; j < rep; j++ {
			var start graph.NI
			for {
				start = graph.NI(r.Intn(len(pos)))
				if len(g.LabeledAdjacencyList[start]) > 0 {
					break
				}
			}
			f := func(b *testing.B) {
				for i := 0; i < b.N; i++ {
					d.AllPaths(start)
					d.Reset()
				}
			}
			b := testing.Benchmark(f)
			fmt.Printf("n=%4d, run %d: %v\n", n, j, b)
			fn := float64(n)
			ln := math.Log(fn)
			l.AddDataPoint(reg.DataPoint{
				Observed: float64(b.NsPerOp()),
				//Variables: []float64{fn, ln, ln * ln, fn * ln, fn * fn},
				Variables: []float64{fn * ln},
			})
		}
		l.RunLinearRegression()
		fmt.Printf("%-7s  %g\n", "C", l.GetRegCoeff(0))
		for i := 0; i < rep; i++ {
			fmt.Printf("%-7s  %g\n", l.GetVarName(i), l.GetRegCoeff(i+1))
		}
		fmt.Println("R2:", l.Rsquared)
	}
	// plot?
	p, err := plot.New()
	if err != nil {
		log.Fatal(err)
	}
	p.Title.Text = "Dijkstra.AllPaths, Directed Graph\nArc Size M = 10N"
	p.X.Label.Text = "Graph Order, N"
	p.Y.Label.Text = "µsec"
	pts := make([]plotter.XYer, len(ns))
	for i, n := range ns {
		xys := make(plotter.XYs, rep)
		pts[i] = xys
		for j := range xys {
			xys[j].X = float64(n)
			xys[j].Y = l.Data[i*rep+j].Observed / 1000
		}
	}
	mmm, err := plotutil.NewErrorPoints(meanMinMax, pts...)
	if err != nil {
		log.Fatal(err)
	}
	if err := plotutil.AddLinePoints(p, mmm); err != nil {
		log.Fatal(err)
	}
	if err := plotutil.AddYErrorBars(p, mmm); err != nil {
		log.Fatal(err)
	}
	p.X.Scale = plot.LogScale{}
	p.Y.Scale = plot.LogScale{}
	p.X.Tick.Marker = plot.LogTicks{}
	p.Y.Tick.Marker = plot.LogTicks{}
	if err := p.Save(4*vg.Inch, 4*vg.Inch, "bench.svg"); err != nil {
		log.Fatal(err)
	}
	// nnls
	A := make([][]float64, len(ns))
	b := make([]float64, len(ns))
	for i, n := range ns {
		fn := float64(n)
		ln := math.Log(fn)
		A[i] = []float64{fn * fn, fn * ln * ln, fn * ln, fn, ln * ln, ln, 1}
		b[i] = mmm.XYs[i].Y
	}
	fmt.Println(nnls.NNLS(A, b, 1000))
}

func meanMinMax(vs []float64) (mean, lowerr, higherr float64) {
	low := math.Inf(1)
	high := math.Inf(-1)
	for _, v := range vs {
		mean += v
		low = math.Min(v, low)
		high = math.Max(v, high)
	}
	mean /= float64(len(vs))
	return mean, mean - low, high - mean
}