Mostly working. Need to check results of tests a bit more.

ackwork/ackwork.go

+// Aug 2021
+
+package ackwork
+
+import (
+	"ackley_mc/ackwork/ackley"
+	"flag"
+	"fmt"
+	"io"
+	"os"
+)
+
+const ndimDflt = 2
+
+type mcPrm struct {
+	iniTmp, fnlTmp float32   // Initial and final temperatures
+	xIni           []float32 // initial x slice
+	xDlta          float32   // change x by up to +- xDlta in trial moves
+	nStep          uint32    // number of steps
+	nRun           uint32    // number of separate runs
+	nOutput        uint32    // write at least this many lines of output
+	fOutName       string    // where we write output to
+	verbose        bool
+	dummy          bool // for testing. If set, do not do a run
+}
+
+var seed int64 // for random numbers, but each thread gets its own value
+
+// usage
+func usage() {
+	u := `[options] input_parameter_file`
+	fmt.Fprintf(flag.CommandLine.Output(), "usage of %s: %s\n", os.Args[0], u)
+}
+
+// realmain is the real main function. The program wants to read from a file,
+// but we can have it read from an io.Reader and then use the test functions
+// for playing with it.
+func realmain(fp io.Reader) error {
+	var mcPrm mcPrm
+	if err := rdPrm(fp, &mcPrm); err != nil {
+		return err
+	}
+ 	if err := doRun(&mcPrm); err != nil {
+		return err
+	}
+	var x [1]float32
+	ackley.Ackley(x[:])
+	return nil
+}
+
+// MyMain is called from the real main. It opens the input file, but passes a
+// reader to the function that does the work. This lets us call it from a test
+// function.
+func MyMain() error {
+	if len(os.Args) < 2 {
+		usage()
+		return fmt.Errorf("No input file given")
+	}
+	fname := os.Args[1]
+	if fp, err := os.Open(fname); err != nil {
+		return err
+	} else {
+		defer fp.Close()
+		if err := realmain(fp); err != nil {
+			return err
+		}
+	}
+	return nil
+}

ackwork/ackwork_test.go

+// Aug 2021
+
+package ackwork
+
+import (
+	"strings"
+	"testing"
+)
+
+func TestErr(t *testing.T) {
+	rdrErr := strings.NewReader("\nrubbish")
+	if err := realmain(rdrErr); err == nil {
+		t.Fatal("line with rubbish should provoke error")
+	}
+}
+
+var s1d = `
+ini_temp 1
+final_temp 0.05
+x_ini 15
+n_step 10000
+n_output 5000
+f_outName testanneal1d`
+
+func TestAnneal1d(t *testing.T) {
+	if err := realmain(strings.NewReader(s1d)); err != nil {
+		t.Fatal("anneal run", err)
+	}
+
+}
+var s3d = `
+ini_temp 1
+final_temp 0.05
+x_ini 15,15,15
+n_step 10000
+n_output 500
+f_outName testanneal3d`
+
+func TestAnneal3d(t *testing.T) {
+	if err := realmain(strings.NewReader(s3d)); err != nil {
+		t.Fatal("anneal run", err)
+	}
+
+}
+
+var sCold = `
+ini_temp 1e-10
+final_temp 1e-12
+x_ini 1,1,1
+n_step 10000
+n_output 500
+f_outName testcold`
+
+func TestCold(t *testing.T) {
+	if err := realmain(strings.NewReader(sCold)); err != nil {
+		t.Fatal("anneal run", err)
+	}
+}

ackwork/dorun.go

+// Aug 2021
+package ackwork
+
+import (
+	"fmt"
+	"math"
+	"math/rand"
+	"io"
+	"os"
+	"sync"
+
+	"ackley_mc/ackwork/ackley"
+)
+
+var seedLocker sync.Mutex
+
+// getSeed returns the random number seed, but we have to put a lock
+// around it so multiple runs get consistent successive seeds
+func getSeed() int64 {
+	seedLocker.Lock()
+	r := seed
+	seed++
+	seedLocker.Unlock()
+	return r
+}
+
+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 ?
+	nEvery   uint32  // Output every nEvery steps
+	fOut     io.WriteCloser
+}
+
+// setupRun does things like get the cooling rate
+func setupRun(mcPrm *mcPrm, cprm *cprm) error {
+	var err error
+	if mcPrm.dummy == true {
+		return nil
+	}
+	cprm.rand = rand.New(rand.NewSource(getSeed()))
+
+
+	if mcPrm.iniTmp != mcPrm.fnlTmp {
+		var coolrate float64
+		cprm.coolme = true
+		if mcPrm.fnlTmp > mcPrm.iniTmp {
+			// would you like to print out a warning ?
+		}
+		nStep := float64(mcPrm.nStep)
+		fnlTmp, iniTmp := float64(mcPrm.fnlTmp), float64(mcPrm.iniTmp)
+		coolrate = -1 / nStep
+		coolrate *= math.Log(fnlTmp / iniTmp)
+		cprm.coolMult = math.Exp(-coolrate)
+	}
+
+	cprm.nEvery = uint32(math.Ceil(float64(mcPrm.nStep) / float64(mcPrm.nOutput)))
+	// change the next line to handle the case of multiple runs
+	ntmp := mcPrm.fOutName + ".csv"
+	if cprm.fOut, err = os.Create (ntmp); err != nil {
+		return err
+	}
+		
+	return nil
+}
+
+// newx gets a candidate x slice
+func newx(xold []float32, xT []float32, rand *rand.Rand, xDlta float32) {
+	for i, x := range xold {
+		t := 2. * rand.Float32() - 1
+		t *= xDlta
+		xT[i] = x + t
+	}
+}
+
+func doRun(mcPrm *mcPrm) error {
+	var cprm cprm
+	if err := setupRun(mcPrm, &cprm); err != nil {
+		return err
+	}
+	if cprm.fOut != nil {
+		defer cprm.fOut.Close()
+	}
+	var nAcc int // number of accepted moves
+	x := make([]float32, len(mcPrm.xIni))
+	xT := make([]float32, len(mcPrm.xIni))
+	nout := uint32(1)
+
+	copy(x, mcPrm.xIni)
+	fOld := ackley.Ackley(x)
+	tmprtr := float64(mcPrm.iniTmp)
+	var n uint32
+	for ; n < mcPrm.nStep; n++ {
+		var acc bool
+		nout--
+		if nout == 0 {
+			fmt.Fprintf (cprm.fOut, "%d,%.4g,%.5g", n, tmprtr,fOld)
+			for _,xx := range x {
+				fmt.Fprintf(cprm.fOut, ",%.2g", xx)}
+			fmt.Fprintln(cprm.fOut)
+			nout = cprm.nEvery
+		}
+
+		newx(x, xT, cprm.rand, mcPrm.xDlta)
+		fTrial := ackley.Ackley(xT)
+
+		if fOld < fTrial {
+			r := cprm.rand.Float64()
+			delta := fTrial - fOld
+			t := math.Exp(-delta / tmprtr)
+			if r < t {
+				acc = true
+			}
+		} else {
+			acc = true
+		}
+		if acc == true {
+			nAcc++
+			copy(x, xT)
+			fOld = fTrial
+		}
+		if cprm.coolme {
+			tmprtr *= cprm.coolMult
+		}
+	}
+	return nil
+}

ackwork/prm_file.go

+// Aug 2021
+// Read the parameters for a calculation
+// We have magic keywords like ini_temp.
+// These are stored in a map along with string values. These are then converted
+// to floats or ints and put into a structure or variable.
+
+package ackwork
+
+import (
+	"bufio"
+	"errors"
+	"fmt"
+	"io"
+
+	"strconv"
+	"strings"
+)
+
+var i int = 6
+var cmdDflt = []struct {
+	key string
+	v   string
+}{
+	{"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", ""},
+}
+
+// digest digests the map, fills out our structure and sets the seed.
+func digest(prmMap map[string]string, mcPrm *mcPrm) error {
+	var err error
+	getf := func(s string) float32 {
+		if err != nil {
+			return 0.0
+		}
+		var r float64
+		r, err = strconv.ParseFloat(s, 32)
+		return float32(r)
+	}
+	getu := func(s string) uint32 {
+		if err != nil {
+			return 0
+		}
+		var u uint64
+		u, err = strconv.ParseUint(s, 10, 32)
+		return uint32(u)
+	}
+	getx := func(s string) []float32 {
+		y := strings.Split(s, ",")
+		if len(y) < 1 {
+			err = fmt.Errorf("broke trying to get initial x values")
+			return nil
+		}
+
+		x := make([]float32, len(y))
+		for i, s := range y {
+			r := 0.
+			if r, err = strconv.ParseFloat(s, 32); err != nil {
+				return nil
+			}
+			x[i] = float32(r)
+		}
+		return x
+	}
+	mcPrm.iniTmp = getf(prmMap["ini_temp"])
+	mcPrm.fnlTmp = getf(prmMap["final_temp"])
+	mcPrm.xDlta = getf(prmMap["x_delta"])
+	mcPrm.nStep = getu(prmMap["n_step"])
+	mcPrm.nRun = getu(prmMap["n_run"])
+	mcPrm.nOutput = getu(prmMap["n_output"])
+	mcPrm.xIni = getx(prmMap["x_ini"])
+	mcPrm.fOutName = prmMap["f_outname"]
+	seed = int64(getu(prmMap["seed"]))
+	if err != nil { // Only returns the first error encountered
+		return err
+	}
+	if prmMap["verbose"] != "" {
+		mcPrm.verbose = true
+	}
+	if prmMap["dummy"] != "" {
+		mcPrm.dummy = true
+	}
+	return nil
+}
+
+// spew spews out an error and known input parameters
+func spew() string {
+	s := `Input file should be of
+key value
+pairs. Values for x_ini are separated by commas without spaces.
+`
+	for _, ss := range cmdDflt {
+		s += fmt.Sprintf("%s %s\n", ss.key, ss.v)
+	}
+	return s
+}
+
+// rdPrmInner is in a function by itself, to make testing easier
+func rdPrm(fp io.Reader, mcPrm *mcPrm) error {
+	prmMap := make(map[string]string)
+	for _, s := range cmdDflt {
+		prmMap[s.key] = s.v // Put default values into map
+	}
+	scn := bufio.NewScanner(fp)
+	for scn.Scan() {
+		s := scn.Text()
+		if len(s) == 0 {
+			continue // blank lines
+		}
+		s = strings.ToLower(s)
+		ss := strings.Fields(s)
+		if len(ss) < 2 {
+			return fmt.Errorf("Only one field in line: %s\n%s", s, spew())
+		}
+		key := ss[0]
+		if _, ok := prmMap[key]; !ok {
+			return fmt.Errorf("unknown keyword %s\n%s", key, spew())
+		} else {
+			prmMap[key] = ss[1]
+		}
+	}
+	if err := scn.Err(); err != nil {
+		return errors.New(spew())
+	}
+	return digest(prmMap, mcPrm)
+}

go.mod

+module ackley_mc
+
+go 1.16