From f3190b0f55a8a8743ca04141e811550f60880f62 Mon Sep 17 00:00:00 2001
From: "Andrew E. Torda" <torda@zbh.uni-hamburg.de>
Date: Sat, 1 Jan 2022 16:29:42 +0100
Subject: [PATCH] Plotting of function value works well enough. Should add
 plotting of temperature and x-values.

---
 ackwork/ackwork.go      |   1 +
 ackwork/ackwork_test.go |  48 ++++++++++++++--
 ackwork/dorun.go        |  37 ++++++++++++-
 ackwork/export_test.go  |   3 +
 ackwork/plot.go         | 120 ++++++++++++++++++++++++++++++++++++++++
 ackwork/rdprm.go        |   2 +
 doc.go                  |  21 +++++++
 go.mod                  |   5 ++
 go.sum                  |   9 +++
 9 files changed, 239 insertions(+), 7 deletions(-)
 create mode 100644 ackwork/export_test.go
 create mode 100644 ackwork/plot.go
 create mode 100644 doc.go
 create mode 100644 go.sum

diff --git a/ackwork/ackwork.go b/ackwork/ackwork.go
index ffc153e..07ef3a5 100644
--- a/ackwork/ackwork.go
+++ b/ackwork/ackwork.go
@@ -20,6 +20,7 @@ type mcPrm struct {
 	nRun           uint32    // number of separate runs
 	nOutput        uint32    // write at least this many lines of output
 	fOutName       string    // where we write output to
+	pOutName       string    // where we plot to
 	verbose        bool
 	dummy          bool // for testing. If set, do not do a run
 }
diff --git a/ackwork/ackwork_test.go b/ackwork/ackwork_test.go
index 30e1c58..55bd69d 100644
--- a/ackwork/ackwork_test.go
+++ b/ackwork/ackwork_test.go
@@ -58,6 +58,7 @@ ini_temp 0.1
 final_temp 0.02
 x_ini 7,7,7,7,7,7,7,7,7,7
 n_step 10000000
+n_step 10
 n_output 10000
 f_outName testanneal10d`
 
@@ -105,13 +106,52 @@ n_step 100000
 n_output 50000
 f_outName testphase`
 
-var plottable_test = []string{
+var csv_test = []string{
 	s1d, s3d, sHot, sCold, sRealAnneal, sPhaseTrans, s5d, s10d}
 
-func TestPlottable(t *testing.T) {
-	for _, s := range plottable_test {
+func TestCsv(t *testing.T) {
+	for _, s := range csv_test {
 		if err := realmain(strings.NewReader(s)); err != nil {
-			t.Fatal("plottable failed with\n", err, "\nInput", s)
+			t.Fatal("csv failed with\n", err, "\nInput", s)
+		}
+	}
+}
+
+var s1dplot = `
+ini_temp 1
+final_temp 0.05
+x_ini 15
+n_step 10000
+n_output 5000
+f_outName testanneal1d
+plot_outName testplot.svg`
+
+var plotTest = []string{
+	s1dplot,
+}
+
+func TestPlot(t *testing.T) {
+	for _, s := range plotTest {
+		if err := realmain(strings.NewReader(s)); err != nil {
+			t.Fatal("plot test failed with\n", err, "\nInput", s)
+		}
+	}
+}
+
+func TestMakepngName(t *testing.T) {
+	var tdata = []struct {
+		in  string
+		out string
+	}{
+		{"boo", "boo.png"},
+		{"boo.", "boo.png"},
+		{"a.png", "a.png"},
+		{"/boo/foo/a.PNG", "/boo/foo/a.PNG"},
+		{"a/b/c/d.svg", "a/b/c/d.png"},
+	}
+	for _, s := range tdata {
+		if tmp := MakepngName (s.in); tmp != s.out {
+			t.Fatal ("Wanted", s.out, "got", tmp)
 		}
 	}
 }
diff --git a/ackwork/dorun.go b/ackwork/dorun.go
index 6551354..ad23767 100644
--- a/ackwork/dorun.go
+++ b/ackwork/dorun.go
@@ -28,8 +28,12 @@ type cprm struct { // parameters calculated from input
 	rand     *rand.Rand
 	coolMult float64 // multiplier for temperature
 	coolme   bool    // Are we cooling or just doing constant temperature ?
+	plotme   bool    // Are we making output for plotting
 	nEvery   uint32  // Output every nEvery steps
 	fOut     io.WriteCloser
+	plotOut  io.WriteCloser
+	plotx    []float64 // Why float 64 ? Because the plotting libraries
+	ploty    []float64 // generally want this
 }
 
 // setupRun does things like get the cooling rate, seed the random numbers.
@@ -44,7 +48,7 @@ func setupRun(mcPrm *mcPrm, cprm *cprm) error {
 		var coolrate float64
 		cprm.coolme = true
 		if mcPrm.fnlTmp > mcPrm.iniTmp {
-			return fmt.Errorf ("Final temp %g higher than initial temp %g", mcPrm.fnlTmp, mcPrm.iniTmp)
+			return fmt.Errorf("Final temp %g higher than initial temp %g", mcPrm.fnlTmp, mcPrm.iniTmp)
 		}
 		nStep := float64(mcPrm.nStep)
 		fnlTmp, iniTmp := float64(mcPrm.fnlTmp), float64(mcPrm.iniTmp)
@@ -62,7 +66,16 @@ func setupRun(mcPrm *mcPrm, cprm *cprm) error {
 	if cprm.fOut, err = os.Create(ntmp); err != nil {
 		return err
 	}
-
+	if mcPrm.pOutName != "" {
+		mcPrm.pOutName = makepngName(mcPrm.pOutName)
+		if cprm.plotOut, err = os.Create(mcPrm.pOutName); err != nil {
+			return err
+		}
+		cprm.plotme = true
+		n_alloc := mcPrm.nStep / 5 // About a fifth of the number of steps
+		cprm.plotx = make([]float64, 0, n_alloc)
+		cprm.ploty = make([]float64, 0, n_alloc)
+	}
 	return nil
 }
 
@@ -76,6 +89,18 @@ func newx(xold []float32, xT []float32, rand *rand.Rand, xDlta float32) {
 	}
 }
 func breaker() {}
+
+// plotPnt adds a point to the data to be plotted. At the start, we check
+// if we have to plot at all. This is done here, to make the main loop
+// below a  bit shorter.
+func plotPnt(n uint32, fTrial float64, plotx, ploty *[]float64, plotme bool) {
+	if !plotme {
+		return
+	}
+	*plotx = append(*plotx, float64(n))
+	*ploty = append(*ploty, float64(fTrial))
+}
+
 // doRun does a Monte Carlo run. Although single precision is fine for the
 // coordinates and function, we use double precision for the temperature.
 func doRun(mcPrm *mcPrm) error {
@@ -86,6 +111,9 @@ func doRun(mcPrm *mcPrm) error {
 	if cprm.fOut != nil {
 		defer cprm.fOut.Close()
 	}
+	if cprm.plotOut != nil {
+		defer cprm.plotOut.Close()
+	}
 	const accRateIdeal = 0.05              // Arbitrary target of 5% acceptance
 	const nRunAccAdj = 200                 // every n steps, check acceptance rate
 	const maxxDlta = 2.                    // max step size
@@ -100,6 +128,7 @@ func doRun(mcPrm *mcPrm) error {
 	nRunAcc := nRunAccAdj // Every nRunAccAdj, try adjusting the step size.
 	const runMult = 0.99
 	xDlta := mcPrm.xDlta // Adaptable step size
+	plotPnt(0, fOld, &cprm.plotx, &cprm.ploty, cprm.plotme)
 	for n := uint32(0); n < mcPrm.nStep; n++ {
 		var acc bool
 		nout--
@@ -145,6 +174,7 @@ func doRun(mcPrm *mcPrm) error {
 			runAcc = runMult*runAcc + (1.0 - runMult)
 			copy(x, xT)
 			fOld = fTrial
+			plotPnt(n+1, fTrial, &cprm.plotx, &cprm.ploty, cprm.plotme)
 		} else { // update the running estimate of acceptance
 			runAcc = runMult * runAcc
 		}
@@ -152,6 +182,7 @@ func doRun(mcPrm *mcPrm) error {
 			tmprtr *= cprm.coolMult
 		}
 	}
+	err := plotWrite(cprm.plotx, cprm.ploty, cprm.plotOut, cprm.plotme)
 	fmt.Println("n accepted:", nAcc, "of", mcPrm.nStep+1)
-	return nil
+	return err
 }
diff --git a/ackwork/export_test.go b/ackwork/export_test.go
new file mode 100644
index 0000000..e8307e6
--- /dev/null
+++ b/ackwork/export_test.go
@@ -0,0 +1,3 @@
+package ackwork
+
+var MakepngName = makepngName
diff --git a/ackwork/plot.go b/ackwork/plot.go
new file mode 100644
index 0000000..369a135
--- /dev/null
+++ b/ackwork/plot.go
@@ -0,0 +1,120 @@
+// 31 Dec 2021
+// This is for doing some plotting. For the moment, this is in the ackwork package,
+// but maybe it should be in its own. It definitely goes in its own file, since
+// it is tied to one of the chart packages I have tried out.
+
+package ackwork
+
+import (
+	"fmt"
+	"io"
+	"math"
+	"path/filepath"
+
+	"github.com/wcharczuk/go-chart/v2"
+	"github.com/wcharczuk/go-chart/v2/drawing"
+	"gitlab.rrz.uni-hamburg.de/Bae5157/axticks"
+)
+
+// makepngName takes a name. If it does not end in png, append .png to the name.
+// If it ends in something else, replace the ending with png.
+// We have stalinistically decided that we will only write png output.
+// svg is nicer, but the files can become silly big.
+func makepngName(fname string) string {
+	const dotpng = ".png"
+	ext := filepath.Ext(fname)
+	switch ext {
+	case "":
+		return fname + dotpng
+	case ".png":
+		fallthrough
+	case ".PNG":
+		return fname
+	default:
+		return fname[0:len(fname)-len(ext)] + dotpng
+	}
+	return "hello"
+}
+
+// maketicks gives reasonable default tick locations
+func maketicks(axisDscrpt axticks.AxisDscrpt, prcsn int) []chart.Tick {
+	xmin, delta, prcsn := axisDscrpt.Xmin, axisDscrpt.Delta, axisDscrpt.Prcsn
+	rnge := axisDscrpt.Xmax - axisDscrpt.Xmin
+	ntick := int(math.Round((rnge / delta) + 1))
+	t := make([]chart.Tick, ntick)
+	const fstr = "%.*f"
+	t[0].Value, t[0].Label = xmin, fmt.Sprintf(fstr, prcsn, xmin)
+	for i := 1; i < ntick; i++ {
+		t[i].Value = t[i-1].Value + delta
+		t[i].Label = fmt.Sprintf("%.*f", prcsn, t[i].Value)
+	}
+	return t
+}
+
+// plotWrt writes out a plot to the given io.writercloser and
+// closes the file when finished.
+func plotWrite(xdata, ydata []float64, plotOut io.WriteCloser, plotme bool) error {
+	if !plotme {
+		return nil
+	}
+	var xAxisDscrpt, yAxisDscrpt axticks.AxisDscrpt
+	var err error
+	if yAxisDscrpt, err = axticks.Tickpos(ydata); err != nil {
+		return fmt.Errorf("plot making y axis: %w", err)
+	}
+	if xAxisDscrpt, err = axticks.Tickpos(xdata); err != nil {
+		return fmt.Errorf("plot making x axis: %w", err)
+	}
+
+	xaxis := chart.XAxis{
+		Name:  "step",
+		Range: &chart.ContinuousRange{
+			Min: xAxisDscrpt.Xmin,
+			Max: xAxisDscrpt.Xmax,
+		},
+		Ticks: maketicks(xAxisDscrpt, yAxisDscrpt.Prcsn),
+	}
+	yaxis := chart.YAxis{
+		Name: "cost (arb units)",
+		NameStyle: chart.Style{
+			TextRotationDegrees: 360}, // Zero does not work
+		TickStyle: chart.Style{
+			TextRotationDegrees: 0,
+		},
+		AxisType: chart.YAxisSecondary,
+		Range: &chart.ContinuousRange{
+			Min: yAxisDscrpt.Xmin,
+			Max: yAxisDscrpt.Xmax,
+		},
+		Ticks: maketicks(yAxisDscrpt, yAxisDscrpt.Prcsn),
+	}
+
+	graph := chart.Chart{
+		Width: 800,
+		Series: []chart.Series{
+			chart.ContinuousSeries{
+				XValues: xdata,
+				YValues: ydata,
+				Style: chart.Style{
+					StrokeWidth: 5,
+					DotWidth:    0,
+				},
+			},
+		},
+		YAxis: yaxis,
+		XAxis: xaxis,
+
+		Background: chart.Style{
+			Padding: chart.Box{
+				Left:  75,
+				Right: 25,
+			},
+			FillColor: drawing.ColorTransparent,
+		},
+	}
+	if err := graph.Render(chart.PNG, plotOut); err != nil {
+		return fmt.Errorf("While plotting %w", err)
+	}
+
+	return nil
+}
diff --git a/ackwork/rdprm.go b/ackwork/rdprm.go
index bfbb810..6b95135 100644
--- a/ackwork/rdprm.go
+++ b/ackwork/rdprm.go
@@ -30,6 +30,7 @@ var cmdDflt = []struct {
 	{"x_delta", "0.5"},
 	{"seed", "1637"},
 	{"f_outname", "mcrun"},
+	{"plot_outname", ""}, // empty means no plots
 	{"verbose", ""},
 	{"dummy", ""},
 }
@@ -81,6 +82,7 @@ func digest(prmMap map[string]string, mcPrm *mcPrm) error {
 	mcPrm.nOutput = getu(prmMap["n_output"])
 	mcPrm.xIni = getx(prmMap["x_ini"])
 	mcPrm.fOutName = prmMap["f_outname"]
+	mcPrm.pOutName = prmMap["plot_outname"]
 	seed = int64(getu(prmMap["seed"]))
 	if err != nil { // Only returns the first error encountered
 		return err
diff --git a/doc.go b/doc.go
new file mode 100644
index 0000000..fe74335
--- /dev/null
+++ b/doc.go
@@ -0,0 +1,21 @@
+// Dec 2021
+
+// ackley is for teaching. We use Monte Carlo to optimise Ackley's function
+// in N dimensions.
+// To run it
+//   Ackley input_file
+// where input_file is a set of name-value pairs. These are defined with defaults
+// in rdprm.go. At the moment, we have
+/* 
+	{"ini_temp", "20"},
+	{"final_temp", "1"},
+	{"n_step", "1000"},
+	{"n_run", "1"},
+	{"n_output", "500"},
+	{"x_ini", "3,4,5"},
+	{"x_delta", "0.5"},
+	{"seed", "1637"},
+	{"f_outname", "mcrun"},
+	{"verbose", ""},
+	{"dummy", ""},
+*/
diff --git a/go.mod b/go.mod
index 0e218ff..b3e69c3 100644
--- a/go.mod
+++ b/go.mod
@@ -1,3 +1,8 @@
 module ackley_mc
 
 go 1.16
+
+require (
+	github.com/wcharczuk/go-chart/v2 v2.1.0
+	gitlab.rrz.uni-hamburg.de/Bae5157/axticks v0.0.0-20211227145643-dc5ef95d7dad
+)
diff --git a/go.sum b/go.sum
new file mode 100644
index 0000000..520af54
--- /dev/null
+++ b/go.sum
@@ -0,0 +1,9 @@
+github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0 h1:DACJavvAHhabrF08vX0COfcOBJRhZ8lUbR+ZWIs0Y5g=
+github.com/golang/freetype v0.0.0-20170609003504-e2365dfdc4a0/go.mod h1:E/TSTwGwJL78qG/PmXZO1EjYhfJinVAhrmmHX6Z8B9k=
+github.com/wcharczuk/go-chart/v2 v2.1.0 h1:tY2slqVQ6bN+yHSnDYwZebLQFkphK4WNrVwnt7CJZ2I=
+github.com/wcharczuk/go-chart/v2 v2.1.0/go.mod h1:yx7MvAVNcP/kN9lKXM/NTce4au4DFN99j6i1OwDclNA=
+gitlab.rrz.uni-hamburg.de/Bae5157/axticks v0.0.0-20211227145643-dc5ef95d7dad h1:IBCLghidSwZff0/3RuvTLDnHgK9SaZHZ/cXWSjimEt0=
+gitlab.rrz.uni-hamburg.de/Bae5157/axticks v0.0.0-20211227145643-dc5ef95d7dad/go.mod h1:5CNHyqeidRypmIVTnQksGqnmP56Oshw1Zv1nlUezrpQ=
+golang.org/x/image v0.0.0-20200927104501-e162460cd6b5 h1:QelT11PB4FXiDEXucrfNckHoFxwt8USGY1ajP1ZF5lM=
+golang.org/x/image v0.0.0-20200927104501-e162460cd6b5/go.mod h1:FeLwcggjj3mMvU+oOTbSwawSJRM1uh48EjtB4UJZlP0=
+golang.org/x/text v0.3.0/go.mod h1:NqM8EUOU14njkJ3fqMW+pc6Ldnwhi/IjpwHt7yyuwOQ=
-- 
GitLab