diff --git a/ackley/ackley.go b/ackwork/ackley/ackley.go similarity index 100% rename from ackley/ackley.go rename to ackwork/ackley/ackley.go diff --git a/ackley/ackley_test.go b/ackwork/ackley/ackley_test.go similarity index 100% rename from ackley/ackley_test.go rename to ackwork/ackley/ackley_test.go diff --git a/ackley/notes b/ackwork/ackley/notes similarity index 100% rename from ackley/notes rename to ackwork/ackley/notes diff --git a/ackwork/ackwork.go b/ackwork/ackwork.go new file mode 100644 index 0000000000000000000000000000000000000000..ffc153e47715d07ae027cfa6dd5de12df459b314 --- /dev/null +++ b/ackwork/ackwork.go @@ -0,0 +1,69 @@ +// 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 +} diff --git a/ackwork/ackwork_test.go b/ackwork/ackwork_test.go new file mode 100644 index 0000000000000000000000000000000000000000..35847b4539f200058286e1bdf11d2847267bcc96 --- /dev/null +++ b/ackwork/ackwork_test.go @@ -0,0 +1,58 @@ +// 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) + } +} diff --git a/ackwork/dorun.go b/ackwork/dorun.go new file mode 100644 index 0000000000000000000000000000000000000000..7463c440d3298cbd529c226f5a3db0c2013358c6 --- /dev/null +++ b/ackwork/dorun.go @@ -0,0 +1,127 @@ +// 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 +} diff --git a/ackwork/prm_file.go b/ackwork/prm_file.go new file mode 100644 index 0000000000000000000000000000000000000000..0f26b1e8a9e784ac1b4964fa02f4ae81537704d1 --- /dev/null +++ b/ackwork/prm_file.go @@ -0,0 +1,134 @@ +// 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) +} diff --git a/function.docx b/function.docx new file mode 100644 index 0000000000000000000000000000000000000000..068019791a47f07899e754c4cfb4861a951d6db4 Binary files /dev/null and b/function.docx differ diff --git a/go.mod b/go.mod new file mode 100644 index 0000000000000000000000000000000000000000..0e218ffc36fd93ea3eb4d037bcebfe0e7182615d --- /dev/null +++ b/go.mod @@ -0,0 +1,3 @@ +module ackley_mc + +go 1.16