From 25da441fcdfbc7a331c36c9860463d541ce522ae Mon Sep 17 00:00:00 2001
From: "Andrew E. Torda" <torda@zbh.uni-hamburg.de>
Date: Fri, 11 Mar 2022 17:50:22 +0100
Subject: [PATCH] Added equilibration steps.

---
 examples/example2  |  4 ++--
 mc_work/dorun.go   | 32 ++++++++++++++++++++++++++++++--
 mc_work/mc_work.go |  3 +++
 mc_work/rdprm.go   |  2 ++
 ui/input_tab.go    | 13 ++++++++-----
 5 files changed, 45 insertions(+), 9 deletions(-)

diff --git a/examples/example2 b/examples/example2
index b163e64..f24c2db 100644
--- a/examples/example2
+++ b/examples/example2
@@ -1,6 +1,6 @@
 # This one is a constant temperature run
-ini_temp 2.0
-final_temp 2.0
+ini_temp 1.1
+final_temp 1.1
 n_step 5000000
 x_ini 4,4
 x_delta 0.5
diff --git a/mc_work/dorun.go b/mc_work/dorun.go
index 718c12b..e767847 100644
--- a/mc_work/dorun.go
+++ b/mc_work/dorun.go
@@ -224,7 +224,7 @@ func nRunAdj(mcPrm *McPrm, nstep int, runAcc float64, xDlta float64) float64 {
 // values, but the interface to the plot and interface packages mostly
 // want double precision, so we have a lot of float64 that shoudl be
 // float32.
-func DoRun(mcPrm *McPrm) (MCresult, error) {
+func singleRun(mcPrm *McPrm) (MCresult, error) {
 	var cprm cprm
 	broken := MCresult{}
 	if err := setupRun(mcPrm, &cprm); err != nil {
@@ -258,7 +258,7 @@ func DoRun(mcPrm *McPrm) (MCresult, error) {
 	const runMult = 0.99
 	xDlta := mcPrm.XDlta // Step size which might be adjusted on the fly
 	saveStep(&cprm, 0, tmprtr, x, fOld)
-
+	
 	for n := 0; n < mcPrm.NStep; n++ {
 		var accept bool
 		if mcPrm.AdaptStep { // Are we trying adaptive steps ?
@@ -326,5 +326,33 @@ func DoRun(mcPrm *McPrm) (MCresult, error) {
 		Fdata: fdata, Xdata: xdata,
 		NStep: mcPrm.NStep, NAcc: nAcc,
 		Bestfval: bestfval, BestX: bestX,
+		LastX: x,
 	}, nil
 }
+
+// DoRun will do a run with equilibration. It will call singleRun twice.
+// On the first run, we make sure there is no cooling and output is
+// not written or goes to the garbage (io.Discard).
+func DoRun(mcPrm *McPrm) (MCresult, error) {
+	var tmpX []float32
+	if mcPrm.NEquil > 0 {
+		tmpX = make ([]float32, len(mcPrm.XIni))
+		mcEquil := *mcPrm
+		mcEquil.FnlTmp = mcEquil.IniTmp
+		mcEquil.NStep = mcEquil.NEquil
+		mcEquil.NEquil = 0
+		mcEquil.fOutName = ""
+		mcEquil.fPltName, mcEquil.xPltName = "", ""
+		if eqResult, err := singleRun(&mcEquil); err != nil {
+			return eqResult, err
+		} else {
+			copy(tmpX, mcPrm.XIni)
+			copy(mcPrm.XIni, eqResult.LastX)
+		}
+	}
+	mcResult, err :=  singleRun (mcPrm)
+	if tmpX != nil {
+		copy (mcPrm.XIni, tmpX)
+	}
+	return mcResult, err
+}
diff --git a/mc_work/mc_work.go b/mc_work/mc_work.go
index 7d79c2e..8c30fe8 100644
--- a/mc_work/mc_work.go
+++ b/mc_work/mc_work.go
@@ -11,6 +11,7 @@ type McPrm struct {
 	XIni           []float32 // initial x slice
 	XDlta          float64   // change x by up to +- xDlta in trial moves
 	NStep          int       // number of steps
+	NEquil         int       // num steps for equilibration
 	NRun           int       // number of separate runs
 	fOutName       string    // where we write output to
 	fPltName       string    // where we plot to (function values)
@@ -29,6 +30,8 @@ type MCresult struct {
 	NStep, NAcc  int       // number of steps and accepted steps
 	Bestfval     float64   // best function value found
 	BestX        []float32 // location of this best value found
+	LastX        []float32 // Last location. Only used after equilibration.
 }
 
+// breaker is for convenience - set breakpoints when debugging.
 func breaker(...interface{}) {}
diff --git a/mc_work/rdprm.go b/mc_work/rdprm.go
index b616e4e..fd4e05f 100644
--- a/mc_work/rdprm.go
+++ b/mc_work/rdprm.go
@@ -22,6 +22,7 @@ var cmdDflt = []struct {
 	{"ini_temp", "20"},
 	{"final_temp", "1"},
 	{"n_step", "1000"},
+	{"n_equil", "0"},
 	{"n_run", "1"},
 	{"x_ini", "3,4,5"},
 	{"x_delta", "0.5"},
@@ -76,6 +77,7 @@ func digest(prmMap map[string]string, mcPrm *McPrm) error {
 	mcPrm.FnlTmp = getf64(prmMap["final_temp"])
 	mcPrm.XDlta = getf64(prmMap["x_delta"])
 	mcPrm.NStep = geti(prmMap["n_step"])
+	mcPrm.NEquil = geti(prmMap["n_equil"])
 	mcPrm.NRun = geti(prmMap["n_run"])
 	mcPrm.XIni = getx(prmMap["x_ini"])
 	mcPrm.fOutName = prmMap["foutname"]
diff --git a/ui/input_tab.go b/ui/input_tab.go
index ff543d5..4b8f97c 100644
--- a/ui/input_tab.go
+++ b/ui/input_tab.go
@@ -114,19 +114,22 @@ func paramForm(genParams genParams) *widget.Form {
 	xIniItem.HintText = xCoordHint
 
 	nStepItem, nStepReload := intItem(&mcPrm.NStep, "N steps", "number of steps in calculation")
+	nEqulItem, nEqulReload := intItem(&mcPrm.NEquil, "N equilibration", "number of steps for equilibration")
 	seedItem, seedReload := intItem(&mcwork.Seed, "seed random numbers", "seed for random number generator")
 
-	adaptbind := binding.BindBool(&mcPrm.AdaptStep)
+	adaptbind := binding.BindBool(&mcPrm.AdaptStep) // checkbox for adaptive stepsize
 	adaptEntry := widget.NewCheckWithData("adaptive step size", adaptbind)
 	adaptItem := widget.NewFormItem("use adaptive step size", adaptEntry)
 	adaptItem.HintText = "Turn on adaptive moves to try to keep acceptance near 5%"
 	adaptReload := func() { _ = adaptbind.Reload() }
+
 	reloadPTab := func() {
-		iniTmpReload()
-		fnlTmpReload()
-		xDltaReload()
+		iniTmpReload() // The values linked to in the the various fields can
+		fnlTmpReload() // change outside of this function, so reloadPTab
+		xDltaReload()  // is called to refresh all the individual elements.
 		xIniEntry.SetText(fslicestrng(mcPrm.XIni))
 		nStepReload()
+		nEqulReload()
 		seedReload()
 		adaptReload()
 	}
@@ -137,7 +140,7 @@ func paramForm(genParams genParams) *widget.Form {
 
 	r := widget.NewForm( // Lays out the items in the form
 		iniTmpItem, fnlTmpItem, xDltaItem, xIniItem,
-		nStepItem, seedItem, adaptItem, rdFileItem)
+		nStepItem, nEqulItem, seedItem, adaptItem, rdFileItem)
 	r.SubmitText = "start calculation"
 	r.OnSubmit = func() { startrun(genParams, r, mcPrm) }
 	reloadPTab() // does this help ? sometimes the form pops up with entries not ready
-- 
GitLab